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ABSTRACT 


The  method  of  characteristics  is  formulated  for  the  computation  of 
the  supersonic  flow  of  an  inviscid,  reacting  gas  over  a  smooth  three-dimensional 
body.  Various  methods  of  constructing  networks  of  bicharacteristic  lines  are 
examined  from  the  point  of  view  of  numerical  stability  and  accuracy.  A  new 
method  of  forming  the  network,  which  consists  of  projecting  forward  along 
streamlines  from  data  points  on  specified  data  planes,  is  found  to  be  most 
easily  adopted  to  the  particular  requirements  of  nonequilibrium  chemistry. 

The  general  method  was  coded  for  the  IBM  7090  computer  and  the 
program  demonstrated  for  the  case  of  an  ideal  gas.  Calculations  were  made 
for  the  flow  about  a  spherical-tip  15“  half-angle  cone  at  10°  angle  of  attack 
and  a  generalized  elliptical  body  at  zero  incidence.  Since  the  program  yields 
the  pressure  distribution  along  specified  streamlines,  it  is  straightforward, 
in  principle,  to  link  it  to  a  finite-rate  chemistry  stream  tube  program  to  treat 
three-dimensional,  nonequilibrium  flows. 
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Section  1 


INTRODUCTION 

1.  I  POSSIBLE  METHODS  OF  SOLUTION 

The  general  equations  governing  the  three-dimensional,  supersonic 
flow  of  a  reacting  gas  are  highly  nonlinear  and  cannot  be  solved  in  closed  form. 
Thus,  without  resorting  to  many  simplifying  approximations,  the  only  alterna¬ 
tive  is  to  solve  the  equations  numerically. 

The  prime  objective  of  the  present  study  was  to  develop  a  numeri¬ 
cal  procedure  for  determining  the  supersonic  flow  field  about  a  general  three- 
dimensional  body  including  nonequilibrium  chemical  effects.  Several  techniques 
are  available  for  doing  this. 

Most  of  these  methods  make  use  of  the  hyperbolicity  of  the  general 
governing  equations  by  formulating  the  problem  as  an  initial  value  problem 
which  can  be  continued  in  either  time  o.  space.  One  method,  however,  which 
has  recently  been  proposed  for  steady  flow  by  Telenin  and  Tinyakov^  treats 
the  problem  as  a  boundary  value  problem  which  must  be  solved  between  two 
boundaries  (i.e.  ,  the  body  and  the  shock)  subject  to  certain  known  boundary 
conditions.  This  scheme,  which  has  not  yet  been  treated  in  detail  in  the 
literature,  is  somewhat  similar  to  the  direct  method  of  Belotserkovsky  for 
two-dimensional  and  axisymmetric  flow  fields  in  that  polynomial  approximations 
are  used  to  reduce  the  partial  differential  equations  to  a  set  of  ordinary 
differential  equations  which  can  be  integrated  numerically.  Both  subsonic 
and  supersonic  flows  can  be  solved.  However,  since  the  method  can  probably 
not  be  easily  extended  to  include  complex  chemistry  models  and  unsteady  flow 
problems,  it  will  not  be  considered  any  further  here. 

Of  the  remaining  possible  methods  only  two  appear  to  be  capable 
of  providing  sufficient  detail  through  the  shock  layer.  The  first  of  these, 
called  the  finite  difference  approach,  consists  of  replacing  the  original  partial 
differential  equations  by  a  system  of  difference  equations  -  formed  through 
direct  substitution  of  finite  differences  for  the  derivates  -  which  can  then  be 
solved  using  a  step-by-step  numerical  scheme.  In  the  second  technique, 
generally  termed  the  method  of  characteristics,  the  equations  of  motion  are 
transformed  to  a  characteristic  coordinate  system  and  the  derivatives  are 
again  approximated  by  finite  difference  equations.  The  resulting  equations  are 
numerically  integrated  step-by-step  along  characteristic  surfaces  throughout 
the  flow-field. 
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The  finite  difference  techniques  can  be  further  subdivided  according 
to  the  form  in  which  governing  partial  differential  equations  are  written  and 
the  differencing  techniques  employed.  In  general,  the  classifications  would 


be  as  follows: 


Standard  (SFD)  -  If  the  equations  of  motion  are  written 
in  either  the  Lagrangian  or  Eulerian  form  and  the  finite 
difference  approximations  are  substituted  directly,  the 
scheme  is  termed  the  standard  finite  difference  technique. 
Discontinuities  are  handled  bv  imposing  appropriate  jump 
conditions  (e.  g.  ,  Rankine-Hugoniot  conditions  at  shocks). 

Finite  difference  technique  utilizing  artificial  viscosity 

(PFD)  -  A  small  "pseudoviscosity"  term  is  introduced  into 

the  nonviscous  flow  equations.  The  resultant  effect  is  that 

the  solutions  remain  stable  even  in  regions  of  large  gradients 

2 

(e.  g.  ,  near  shock  waves)  Von  Neumann  and  Richtmyer 

obtained  a  solution  to  the  one -dimensional,  unsteady  flow 

equations  using  this  method  which  accounts  for  the  presence 

of  free  boundaries  (e.  g.  ,  shockwaves)  automatically, 

3 

Burstein  has  recently  extended  this  idea  to  two-dimensional, 
unsteady  flows. 

Finite  difference  procedure  applied  to  conservation  laws 
4 

(CFD)  -  Lax  developed  a  method  which  treats  the  governing 

equations  in  conservation  form  and  uses  a  differencing 

procedure  which  has  the  effect  of  introducing  dissipative 

terms.  The  method  can  handle  cases  in  which  the  solution 

has  certain  specific  kinds  of  discontinuities  (termed  weak 

5  3 

solutions).  Lax  and  Wendroff  ,  Burstein  and  Bohachevsky 
have  recently  applied  this  method  to  various  fluid  flow 
problems.  Also,  Moretti^  has  used  this  procedure  coupled 
with  a  method  of  characteristics  approach  at  the  boundary 
points  to  solve  unsteady  flow  equations. 
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4.  Finite  difference  procedure  applied  directly  to  Navier-Stokes 
equations  (NSFD)  -  This  procedure  which  is  generally 
difficult  because  of  the  complexity  of  the  differential  equations 
involved  was  first  applied  by  Ludford  et  al.  ,  and  more 
recently  by  Crocco  for  one-dimensional  unsteady  flow. 

While  other  methods  have  been  proposed,  these  four  appear  to  be  the  most 
useful  for  calculating  multidimensional  flow  fields. 

In  Table  I,  a  summary  of  some  of  the  advantages  and  disadvantages 

of  each  of  these  finite  difference  approaches  and  the  method  of  characteristics 

7  8  9 

(MOC)  approach  are  given.  ’  ’  On  the  basis  of  these  comparisons  the 

following  conclusions  appear  warranted. 

1.  Although  finite  difference  methods  can  treat  shock  waves 
and  other  discontinuities,  their  accuracy  is  impaired 
(especially  near  the  shock)  because  they  have  the  effect 
of  smearing  out  the  shock  wave  over  several  mesh  points. 

The  effect  of  this  smearing  on  the  rest  of  the  flow  field 
for  multidimensional  flows  is  not  known  at  present. 

2.  Fixed  boundaries  are  more  accurately  treated  by  the 
MOC  technique  when  the  body  shape  is  arbitrary. 

The  finite  difference  schemes  appear  to  be  better  only 
when  the  body  shape  is  such  that  points  on  the  body  lie 
along  coordinate  lines  so  that  body  points  are  mesh 
points.  Otherwise,  interpolations,  which  could  lead  to 
errors  and/or  possibly  instability,  are  required. 

3.  The  disadvantages  of  the  method  of  characteristics  approach 
are  that  the  program  logic  required  to  code  the  procedure 
for  the  computer  is  quite  complicated  and  that  in  order  to 
obtain  a  solution  at  a  given  mesh  point,  considerable  iterating 
is  needed.  This  latter  situation,  of  course,  adversely 
affects  the  computing  time  required  to  obtain  a  solution. 
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TABLE  ! 

COMPARISON  OF  POSSIBLE  APPROACHES 


MOC 

SFD 

PFO 

CFD 

NSFD 

1.  FREE 

BOUNDARIES 

(SHOCKS, 

CONTACT 

SURFACES, 

ETC.) 

TREATED  AS 
DISCONTIN¬ 
UITIES  USING 
SPECIAL 
RELATIONS 

TREATED  AS 
DISCONTIN¬ 
UITIES  USING 
SPECIAL 
RELATIONS 

AUTOMATICALLY 
HANOLED  -  HOW¬ 
EVER  RESOLUTION 

1 S  SPREAO  OVER 
SEVERAL  MESH 
WIDTHS 

SAME  AS 

PFO 

SAME  AS 

PFO 

2.  FIXED 

BOUNDARIES 

EXACT 

BOUNDARY 

IS  SOLVED 

USING 

SPECIAL 

relation 

SAME  AS 

MOC 

DIFFICULTY  CAN 
ARISE  BECAUSE 
ARTIFICIAL 
VISCOSITY  IN¬ 
CREASES  ORDER 

DF  EQUATIONS 

TREATED 

BY  INTER¬ 
POLATION 

TREATED 

BY  SPECIAL 
RELATIONS 

THUS  REQUIRING 
EXTRA  BOUNDARY 
CONDITIONS 


3. 

NO.  OF  INDE¬ 

3-SPACE 

SAME  AS 

SAME  AS  HOC  BUT 

SAME  AS 

SAME  AS 

PENDENT 

PLUS 

MOC 

MORE  0 1 FF 1  CULT 

MOC 

MOC 

VARIABLES 

TIME 

BECAUSE  EXTRA 

WHICH  CAN  BE 

BOUNDARY  CONDI¬ 

SOLVED 

TIONS  MUST  BE 

INTRODUCED 

1. 

STEP  SIZES 

POTENTIALLY 

SMALLER  THAN 

SMALL  BECAUSE 

same  AS 

SAME  AS 

LARGEST 

MOC  BECAUSE 

HIGH 

PFD 

PFD 

TRUNCATION 

GRADIENTS 

ERRORS  ARE 

REQUIRE  A 

LARGER 

FINER  MESH 

5. 

SOLUTION 

UNEQUALLY 

EQUALLY 

SAME  AS 

SAME  AS 

SAME  AS 

GRIO 

SPACED 

SPACED 

SFD 

SFD 

SFD 

6. 

ACCURACY 

POTENTIALLY 

INTERMEDIATE 

LEAST  ACCURATE 

ACCURACY 

SAME  AS 

MOST  ACCU¬ 

ACCURACY 

BECAUSE  ERRORS 

IS  IM¬ 

CFO 

RATE  BECAUSE 

BECAUSE  DF 

CAN  BE  INTRO¬ 

PAIRED  BY 

THIS  METHOD 

LARGER 

DUCED  BY  NOT 

smearing 

OF  SOLUTION 

TRUNCATION 

PROPERLY  SPECI¬ 

OUT  SHOCK 

MOST  CLOSELY 

ERRORS 

FYING  EXTRA 

OVER 

FOLLOWS  FLOW 

BOUNDARY 

SEVERAL 

MODEL 

CONDITIONS 

MESH  POINTS 

7. 

PROGRAMING 

BY  FAR  MOST 

OF  INTER¬ 

SIMPLEST 

SAME  AS 

SAME  AS 

DIFFICULTIES 

DIFFICULT 

MEDIATE 

PFD 

PFD 

DIFFICULTY 
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As  noted  initially,  the  prime  objective  of  the  present  study  was  to 
develop  a  numerical  procedure  for  determining  the  flow  field  about  a  general 
3-D  body  including  nonequilibrium  chemical  effects.  Particular  interest 
was  centered  in  the  afterbody  or  supersonic  portion  of  the  flow  field  because 
for  a  typical  spherically  blunt  nose,  the  forebody  or  subsonic -transonic 
portion  of  the  flow  field  can  be  obtained  from  a  transformation  of  an  appro¬ 
priate  axisymmetric  solution  except  at  high  angles  of  attack.  One  naturally 
wonders  whether  this  problem  can  be  most  efficiently  treated  using  the 
method  of  characteristics  or  one  of  the  finite-difference  approaches.  The 
answer  to  the  question  depends  on  three  considerations: 

(a)  Computing  time  per  solution 

(b)  accuracy  and  reliability  of  the  solution  near  the  shock  and 
body 

(c)  prospects  for  treating  chemical  nonequilibrium  effects. 

When  the  present  work  was  started  in  June  1963,  there  was  not  sufficient 
information  available  concerning  the  finite-difference  methods  to  even  begin 
to  consider  the  answer  to  this  question  rationally.  In  fact  at  the  present 
time,  after  3-4  years  experience  with  both  general  methods,  the  answer  to 
this  question  is  still  not  clear  in  the  case  of  reacting  gas  flows.  In  1963,  it 
was  decided  to  pursue  the  method  of  characteristics  because  considerable 
success  had  been  attained  using  that  method  for  2-D  and  axisymmetric  flows 
and  it  yielded  results  having  the  desired  detail  near  the  shock  and  body.  It 
was  also  clearer  how  procedures  for  treating  nonequilibrium  chemical  effects 
could  be  included  in  the  program  to  various  degrees  of  approximation  (ex. 
chemically  frozen  along  streamlines,  finite-rate  streamtube  with  locally 
frozen  gas  chemistry  for  obtaining  flow  quantities  in  the  mesh  calculation, 
etc.  ).  Hence  this  report  deals  in  detail  with  applying  the  method  of  char¬ 
acteristics  to  three-dimensional  reacting  flows  and  no  further  consideration 
is  given  to  finite-difference  methods.  Questions  regarding  the  relative  merits 
of  the  two  approaches  must  await  further  work. 
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1.2 


SOLUTION  OF  HYPERBOLIC  PARTIAL  DIFFERENTIAL 
EQUATIONS  USING  THE  METHOD  OF  CHARACTERISTICS 


The  theory  of  the  method  of  characteristics  as  applied  to  quasi  - 
linear,  hyperbolic  equations  has  been  known  since  Monge  and  Ampere  in 

g 

the  early  nineteenth  century.  A  monograph  by  Massau  first  indicated  the 

use  of  the  theory  to  solve  a  system  of  two  equations  in  two  unknowns.  Later, 

in  1928,  Lewy^used  the  method  to  show  that  the  initial  value  problem  for 

a  quasi-linear  hyperbolic  differential  equation  in  two  independent  variables 

has  a  unique  solution.  Titt^  later  generalized  the  method  of  Lewy  to  the 

case  of  three  independent  variables.  Since  that  time,  the  method  has  been 

utilized  by  several  authors  in  formulating  methods  for  solving  many  different 

12  13 

hyperbolic  problems  in  fluid  mechanics.  Ferri  and  Meyer  discuss  much 
of  the  earlier  work  as  it  applied  to  several  different  problems  in  fluid 
dynamics  and  give  a  complete  list  of  references.  Now,  with  the  recent 
development  of  faster  and  larger  electronic  digital  computers,  more 
ambitious  problems  such  as  the  one  described  here,  may  be  attempted. 

Thus,  it  is  now  possible  to  consider  the  flow  over  general  bodies  exhibiting 
rather  general  motions  and  including  complicated  thermochemical  effects. 

The  Cauchy  problem  (or  initial  value  problem)  considered  here 
can  be  summarized  simply  as  follows.  Given  an  initial  data  surface  and 
known  boundary  conditions  (body  surfaces  or  shock  waves),  the  method  of 
characteristics  is  to  be  used  to  obtain  the  solution  of  the  governing  partial 
differential  equations  at  another  surface  separated  in  either  time  or  space  from 
the  original  one.  The  solution  will  be  obtained  at  a  discrete  set  of  grid 
points  on  the  new  surface  by  integrating  the  governing  equations  numerically 
along  certain  selected  characteristic  lines.  When  as  many  new  grid  points 
as  needed  are  obtained,  the  surface  thus  calculated  is  taken  as  the  new  initial 
value  surface  and  the  next  adjacent  surface  is  obtained  similarly.  The 
process  is  then  continued  step-by-step  until  the  desired  portion  of  the  flow- 
field  is  solved.  This  technique  has  been  used  successfully  by  innumerable 
authors  for  solving  flow  fields  in  two  independent  variables.  However,  very 
little  actual  numerical  work  has  been  performed  for  problems  involving  either 
complicated  thermochemical  effects  or  more  than  two-independent  variables. 
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It  is  the  object  of  the  present  work  to  obtain  a  numerical  method  for  solving 
just  such  problems. 

Several  authors  have  discussed  the  applications  of  the  method  of 

characteristics  to  multidimensional  flow  fields  from  a  purely  theoretical 

viewpoint  without  offering  any  numerical  solutions.  These  include:  (1)  the 

14 

pioneering  efforts  of  Thornhill  ,  who  discussed  two  possible  characteristic 
networks  for  three-independent  variables,  (2)  Coburn  and  Dolph  ,  who 

1  5 

extended  the  work  of  Titt  to  fluid  flow  problems,  (3)  Clippinger  and  Giese  , 

who  used  a  generalized  vector  formulation  to  derive  the  basic  characteristic 

equations  and  (4)  Holt^,  who  established  a  finite  difference  scheme  based 

17 

upon  the  earlier  work  of  Coburn  and  Dolph.  More  recently,  Fowell  has 
written  the  basic  finite  difference  equations  which  apply  to  one  of  the  networks 
which  Thornhill  proposed.  Although,  as  will  be  discussed  in  Section  3. 
Fowell's  method  was  shown  to  be  unstable,  it  has  proved  to  be  very  useful 
in  developing  stable  multidimensional  integration  analogs. 

A  second  group  of  authors  have  succeeded  in  obtaining  some  results 

using  numerical  computational  procedures.  The  first  few  of  these  were 

limited  to  simpler  calculations  which  were  accomplished  by  hand  with  the 

1 8 

aid  of  desk  calculators.  Among  these  are  included  the  works  of  Ferrari  , 

19  20  21  22 

Moeckel  and  Bruhn  and  Haack  .  Later  Butler  '  and  Tsung  succeeded 

in  obtaining  some  limited  solutions  by  utilizing  the  earlier  digital  computers. 

23  24  9 

More  recently,  Moretti  et  al.  ,  Kackova  and  Cuskin  and  Sauerwein 

have  published  results  obtained  with  the  present  generation  of  high  speed 

computers.  The  most  important  of  these  various  schemes  will  be  discussed 

and  compared  in  Section  3.  It  would  appear  that  the  possibilities  in  this 

area  have  just  begun  to  be  fully  realized  and  with  the  addition  of  still  faster 

and  more  efficient  computers  even  more  complicated  problems  should  prove 

amenable  to  solution  using  the  method  of  characteristics  approach. 

The  extensions  required  to  the  MOC  to  include  nonequilibrium 

Z  5  Z  6 

chemistry  effects  have  been  set  forth  by  Chu  ,  Resler  and  Wood  and 
2?  29 

Kirkwood  among  others.  Sedney  et  al.  have  applied  the  numerical  method 

to  calculate  vibrational  relaxation  and  Capiaux  and  Washington"^  and 
3  1 

Eastman  have  utilized  a  simDlified  reaction  system  without  relaxation. 
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Zupnik  et  al.  ,  and  Widawsky  have  recently  developed  and  applied  the 

34 

theory  for  calculating  nonequilibrium  nozzle  flows,  while  only  Wood  et  al. 

35 

and  Curtis  have  applied  the  method  to  obtain  numerical  results  for  the 
afterbody  flow  fields  including  the  full  reacting  and  relaxing  flow  equations. 

The  present  effort  describes  the  development  of  a  method  which 
in  principle  is  capable  of  calculating  the  flow  field  in  the  vicinity  of  a  general 
three-dimensional  body  including  nonequilibrium  thermochemical  effects. 
Since  the  development  includes  the  unsteady  flow  equations,  the  scheme  could 
eventually  be  used  to  calculate  the  flow  properties  in  the  vicinity  of  a  blunt 
asymmetric  body  or  a  sharp  cone  at  angle-of-attack.  In  Section  2  the  basic 
theory  of  the  method  of  characteristics  as  it  applies  to  the  present  problem 
is  given.  The  relations  are  written  in  terms  of  two  special  characteristic 
coordinate  systems,  each  of  which  is  thought  to  be  best  for  numerical 
applications.  The  proposed  integration  network  is  presented  in  Section  3 
and  is  discussed  and  compared  with  other  schemes  which  have  previously 
been  utilized.  The  discussion  includes  the  requirements  which  are  thought 
necessary  for  insuring  stability  of  a  successful  integration  procedure.  A 
description  of  the  numerical  scheme  as  it  applies  to  each  of  three  unit 
processes  for  calculating  the  flow  properties  at  points  in  the  field,  on  the 
body  surface  and  on  the  shock  wave  surface  is  given  in  some  detail,  including 
the  modifications  necessary  to  account  for  nonequilibrium  chemical  effects. 

The  remainder  of  the  work  describes  the  actual  machine  program 
in  which  the  methods  proposed  in  Section  3  are  applied  to  solve  practical 
flow  problems.  In  the  present  case,  the  program  is  limited  to  steady,  super¬ 
sonic,  three-dimensional,  frozen-inhomogeneous  or  ideal  gas  flows;  however, 
the  program  was  written  in  such  a  way  that  the  extension  to  more  complicated 
thermochemical  models  is  relatively  easily  accomplished.  A  discussion  of 
the  problems  encountered  in  writing  the  program  and  some  typical  results 
for  the  flow  fields  around  conical  and  elliptical  afterbodies  are  given. 

A  description  of  the  machine  program  and  comments  concerning 
the  functions  of  each  subroutine  are  presented  in  Appendix  A.  Flow  diagrams, 
output  description  and  input  and  operating  procedures  are  also  included. 
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Section  2 


DERIVATIONS  OF  GENERAL  EQUATIONS 

In  this  section  the  equations  which  are  required  to  solve  a 
general  nonequilibrium,  three-dimensional  steady  flow  problem  numerically, 
utilizing  the  method  of  characteristics,  are  derived.  The  general  theory 
of  the  method  of  characteristics  as  it  applies  to  a  system  of  first  order 
hyperbolic  partial  differential  equations  is  first  developed.  Then,  the  general 
equations  governing  the  flow  of  a  reacting  and  relaxing  gas  are  presented. 
Finally,  the  theory  will  be  applied  to  the  given  governing  equations  in  order 
to  determine  the  characteristic  conditions. 


2.  1  THE  METHOD  OF  CHARACTERISTICS 

The  general  theory  for  the  method  of  characteristics  is  well 

developed  and  may  be  considered  classical.  Only  that  part  of  the  theory 

which  is  necessary  for  the  development  of  the  proposed  numerical  scheme 

will  be  presented  here.  A  more  complete  development  may  be  found  in 

many  works  on  partial  differential  equations  among  which  the  more  notable 
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are  Courant  and  Hilbert  and  Hadamard 

A  general  system  of  n  -first  order,  quasi-linear  partial 
differential  equations  in  n  -dependent  variables  u and  m  -  independent 
variables  X-  is  considered.  Using  the  index  summation  notation  of  cartesian 
tensors,  these  can  be  written  in  the  form 


da 


(1) 


where  the  <2^^  and  are  given  functions  of  the  and  x>k  .  Quasi- 
linear  systems  of  equations  are  those  in  which  the  highest  order  derivatives 
appear  linearly.  The  restriction  here  to  first  order  equations  does  not 
imply  any  loss  of  generality,  because  systems  containing  higher  order 
derivatives  can  be  reduced  to  a  system  of  first  order  equations  by 
defining  new  dependent  variables. 
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The  method  of  characteristics  in  two-independent  variables 
which  transforms  the  governing  partial  differential  equations  into 
ordinary  differential  equations  along  certain  characteristic  directions  is 
a  very  special  case  of  the  general  theory,  and  hence  it  cannot  be  easily 
generalized.  However,  the  basic  property  that  characteristic  curves  in 
two-variables  are  curves  along  which  the  derivatives  are  continuous  but 
across  which  derivative  discontinuities  may  occur  can  also  be  taken  as 
the  definition  of  a  characteristic  in  m  -independent  variables.  This 
property  is  used  in  two-dimensions  in  order  to  identify  the  characteristics 
as  paths  of  waves  in  the  physical  model.  In  general  then,  a  characteristic 
in  m  -independent  variables  can  be  defined  as  a  subspace  of  m-  1  - 
dimensions  along  which  the  derivatives  of  the  dependent  variables  are 
continuously  differentiable  but  across  which  discontinuities  in  the  variables 
are  allowed  to  occur,  A  characteristic  surface  (or  hypersurface  in  more 
than  three-independent  variables)  is  thus  once  again  associated  with  the 
surface  generated  by  a  wave  front. 

If  the  dependent  variables  are  given  as  smooth  functions  on 
a  characteristic  surface,  it  then  follows  that  the  "interior"  derivatives, 
that  is,  derivatives  tangential  to  a  characteristic  are  continuous  and  only 
the  "exterior"  derivatives  (derivatives  normal  to  the  surface)  may  be 
discontinuous . 

We  consider  a  characteristic  surface  whose  equation  is 


f(xt)  =  0  (2) 

and  apply  a  transformation  to  the  m  -independent  characteristic  variables 
)0t'  by  taking 


V-  H*i)  i  *z  *'}=  *(*i)  ,  V" 


(3) 


Here,  m  -  4  and  the  functions  y,  h  and  k  are  arbitrary  provided  they 

form  an  independent  set  with  f  .  At  the  characteristic  -  0  the 

J  3  a 

derivatives  ✓  and  — -  are  "interior"  to  the  characteristic 


Jx/  >  <?*3' 


and  are  therefore  continuous  (since  u,^  are  smooth).  Thus,  only  the 

derivatives  rnay  be  discontinuous. 

e>x; 
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Now,  the  derivative  of  any  dependent  variable  U p  with  respect  to  one  of  the 

original  independent  variables  is  related  to  its  derivatives  with  respect  to 

/ 

the  new  variables  (  X  •  )  by  the  equations, 


du. p  3  x>j  du>p 
<?zt  ' 

The  jump  in  at  %>/  -  0  is  therefore 

<XXit 

3Up  3oCp  <9 ci.^,  3f 

3xt  3x/ 


(4) 


(5) 


since  all  other  derivatives  are  continuous  as  mentioned  above.  If  (5)  is 
substituted  into  the  set  of  equations  (1)  a  set  of  linear  homogeneous 
equations  for  J  with  coefficients  depending  upon  the  -  and  the 

dependent  variables  Up  is  obtained  as 


<  U  j) 

where  the  u  i  are  transformed  coefficients  of  .  which  are 

known  functions  of  the  Up  and  Xf'  .  Not  all  of  the  jumps  j  can  be 

zero  if  a  genuine  discontinuity  is  possible.  Hence,  the  determinant  of  the 

/ 

coefficient  array  must  v  anish,  giving  a  relationship  between  the 

coefficients  which  must  be  satisfied  if  -f  (x-)~  0  is  to  be  a  characteristic. 
Thus,  the  characteristic  surfaces  (or  hypersurfaces)  are  defined  by 

=  dtt  I °  (7) 

Equation  (7)  is  generally  termed  the  characteristic  equation  for  the 
characteristic  x>‘-  constant. 

The  vanishing  of  the  coefficient  determinant  (7)  can  be  seen  to 
have  another  implication.  Writing  the  original  equations  (1)  in  terms  of 
the  new  independent  variables  (?/,  *t' ,  jsf',  *  *  ')  a  set  of  equations  is 
obtained  which  will  have  the  form 
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(8) 


Where,  the  functional  form 


-  Q'/Li.yi. 


(9) 


is  the  matrix  of  the  system  of  transformed  equations  correspond  ng  to  (1) 
when  /  is  replaced  by  $  ,  and  the  vector  Bp.  contains  all  of  the  terms  in 
the  original  equations  which  do  not  contain  derivatives. 

Now,  if  f  ( )  is  to  be  a  characteristic  manifold  then  by  (7)  the 
determinant  of  Apy  (/  )  must  vanish.  Thus  the  rows  of  (  r  )  must 
be  linearly  dependent  and  there  exists  a  vector  h p  such  that 


*  o  on  (10) 


If  Equations  (8)  are  multiplied  by  this  vector  a  linear  combination  of  the 

equations  is  obtained  which  does  not  involve  derivatives  with  respect  tc 
/ 

X>,  as 


hp.AjU't)  (q)  ^  (h)j£-r-  +  *  -  (H) 

Equation  (11)  is  called  the  characteristic  or  compatibility  condition  and 
does  not  involve  derivatives  with  respect  to  the  coordinate  X>t'  ,  i.e.  ,  in 
the  direction  normal  to  a  characteristic  surface.  There  may  be  more  than 
one  such  equation  associated  with  a  particular  characteristic  surface,  but 
there  must  always  be  at  least  one.  Reference  38  contains  a  discussion  of 
the  possible  number  of  compatibility  equations  associated  with  a  given 
characteristic  manifold.  The  compatibility  equations  are  fundamental  to 
the  method  of  characteristics  and  will  be  derived  below  for  the  cases  of 
three-dimensional  steady  and  unsteady  fluid  flow. 
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2.2 


THE  EQUATIONS  OF  MOTION 


The  equations  of  conservation  of  species,  mass,  momentum 
and  energy  along  with  an  equation  of  state  govern  a  three-dimensional 
nonsteady,  nonequilibrium  gas  flow.  Neglecting  all  transport  effects 
(viscosity,  heat  conduction,  etc.)  these  equations  in  non-dimensional 
form  become: 

Continuity  Equation 


dp  dp  dp  ^  dp  / da  d  ir  dur  \ 
dt  dz  dy  d}  ~\dx  dy  dy  / 


dy  d$ 


(12) 


X  -  Momentum 


da  da  da  da  1  dp 

+  U,  +  ir  —  +  ur  —  +  —  -f— 

di •  dz  dy  dj.  p  dz 


-  O 


(13) 


y  ~  Momentum 


dir  dir  dir  dir  1  dp 

—  +  a  — r  V  —  +  ur - r - — 

dt  dx.  dy  d}  yo  dy 


=  O 


(14) 


-  Momentum 


dur  dur  dur  did  f  dp 

- +  a  —  /  ir  —  +  ur  - —  +  —  — 

dt  dz  dy  dj.  p  d^ 


(15) 


Energy  Equation 


dh  d/>  dd  dh  -A. /dp  dp  dp  dp  1 

dt  dp  dy  df.  p\dt  d%  dy  di' 


dy  d l 


(16) 


Conservation  of  *  th  Species 


i£  ^4* 

dt  dz  dy 


Q 

=  *  -  £+/.  •  •  ,NS  (17) 

/> 
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Conservation  of  i  th  Element 


+  u. 


dx 


*¥ 


+  us 


£4 


0  k-t,2.--,£  (18) 


Vibrational  Energy  Equation 

3e-  ,  Sej  d&. 

— i  +  u  — <L  -h  is  — L+  us  —-i 

tfx>  Jy 


=  *„.  J  -  /,  NV  (19) 

tl 


Equations  of  State 


NS 

h  sZ  \  K 

i  *  / 


and 


ns 

A P  S/>TZ  Yi 

i:/ 


(20) 


(21) 


Equations  (1  2)-(21)  represent  MS*/VV+  7  equations  in  tJS+NV+7  unknowns 
with  NS+NV-£  functions  (  ,  Rv^  )  to  be  specified  and  thus  form  a 

complete  set.  Most  of  the  thermodynamic  and  chemical  variables  have 
been  non-dimensionalized  by  using  their  free-stream  counterparts.  The 
exceptions  to  this  are  the  pressure  which  is  normalized  with  respect  to 
Um  and  the  energies  for  which  is  used.  Distances  are 

normalized  with  respect  to  a  characteristic  length  L?  (which  in  the  case  of 
a  blunted  body  is  usually  the  nose  radius  of  curvature),  and  time  with 

l' 

respect  to  — —  .  In  the  Equation  of  State  (21)^  is  the  normalizing  term 

J  oo 


.A  = 


/te 


(22) 


The  gas  in  the  shock  layer  has  been  assumed  to  consist  of 
M 5  monatomic  and  diatomic  species  Nfj  connected  by  the  R  reactions 
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NS  kf ■ 

./  -  f  j*f 


i-  1,2.-  -.NR  (23) 


where  if- •  and  if- - *  represent  the  stoichiometric  coefficients. 
ij  tj 

Concentrations  are  expressed  in  mole  fractions  referred  to  the  free  stream 
molecular  weight  or 


*  SkL 

i*  i 


Equations  (1  7)and (1  8) govern  the  rate  of  change  of  the  concentration  of  the 
i  th  species.  Here,  QC  ^  is  the  number  of  atoms  of  the  k  th  element 
per  molecule  of  the  i  th  species,  £  represents  the  number  of  elements 
and  NS  represents  the  total  number  of  species.  Equation  (21)  governs  the 
rate  of  vibrational  energy  relaxation  where  £ j  is  the  vibrational  energy 
for  each  of  the  NV  species  which  are  not  in  vibrational  equilibrium. 

The  source  terms  <0-  and  /?„.  are  complicated  functions  of 
the  temperature,  density  and  the  rates  of  certain  reactions  involved  in 
the  species  i  and  j  .  Since  the  exact  forms  of  the  source  terms  are 
not  necessary  for  a  development  of  the  proposed  numerical  solution,  a 
full  discussion  of  them  shall  not  be  given  here.  Reference^  (39)  and  (33) 
describe  the  source  terms  which  are  currently  being  used  for  non¬ 
equilibrium  flow  field  calculations. 

The  enthalpy  of  the  J  th  species  used  in  Equation  (20)  and  the 
other  species  thermodynamic  quantities  are  generally  calculated  from  a 
single  harmonic  oscillator  approximation  to  the  diatomic  molecules  or 
from  curve  fits  as  functions  of  temperature.  Here,  the  caloric  equation 
of  state  of  the  L  th  species  is  given  as 


hL  =  XL  +  +  CL 
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where  includes  the  energy  of  rotation  and  translation,  h L  the  neat  of 
formation  and  6t  the  energy  of  vibration.  Equations  (25),  (20)  and  (21) 
can  be  used  to  eliminate  the  derivatives  of  ft  in  Equation  (12)  Differentiat¬ 
ing  Equation  (20)  gives 


dh 

dt 


i-t 


(26) 


d  3  ,  eP  (?  c? 

where  ~  +  ^  Ifx  ~chf.  +  67 ^  is  "su6stantial"  or 

Stokes  derivative.  Substituting  for  A6  from  (25) 


/  t  NS  /  j  -  \  _/  T"  y  j .  . .  j  -r 

i~  /  4*  t  '  /  ,/■=•' 


Then  by  defining 


ns  ,r  Atf-AV  «- 

£  =  r  y  -  El  +  r  y  El 

Op  /_  rt  dT  4~  dT 

ft  4  *  / 


as  the  specific  heat  at  constant  pressure,  Equation  (27)  becomes 


i±.c  f 

«  '  L  dt  £  h‘  dt  £  *  a 


The  equation  of  state  (21)  can  be  used  to  eliminate  from  Equation  (29) 

Differentiating  (21)  and  combining  the  result  with(29)  gives 


'  Cg*UN  d£  _  L  <*£  .  —  -  -/— )  Y  2i±i  ,  yrdv 

P  c.Xr  dt  i°  dt  L,  dt  (crrJ  P  LY‘“V' 

t  i  ml 


where  fi  =  —  and  hW  = 


NS 

z  *.■ 
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Eliminating  the  term  from  (12)  and  (30)  gives  finally 

'  cLt 


1  dp 
r jO a,*  cU 


7't  +7T‘  m  ■  -f'f 


(32) 


where 


J  /  NV  ns  \  US  ^ 


(33) 


contains  all  the  nonequilibrium  rate  terms  from  Equations  (17)  and  (19)  and 


Cp.  *7 
Cfi  AM/-1 


(34) 


25 

is  the  frozen  speed  of  sound.  Chu  has  shown  that  the  frozen  sound 
speed  a,  defined  by  (34),  is  the  wave  velocity  of  a  flow  which  is  not  in 
equilibrium  and  thus,  can  be  associated  with  the  rate  of  propagation  of 
characteristic  surfaces.  This  result  will  be  derived  in  the  next  section. 
Equations  (13),  (19)  and  (33)  become  the  basic  equations  of  change  to  which 
the  theory  of  characteristics  developed  in  the  last  section  will  be  applied. 


2.  3  THE  CHARACTERISTIC  RELATIONS 

Here,  the  theory  presented  in  Paragraph  2.  2  is  applied  to  the 
flow  equations  given  in  Paragraph  2.  1.  The  governing  equations  for 
steady  flow  are  obtained  from  Equations  (12)  through  (21)  by  setting  the 
terms  involving  ^  identically  equal  to  zero. 

Applying  Equation  (5)  to  the  governing  flow  equations  presented 
in  Paragraph  2.  2  leads  to  the  following  homogeneous  equations  for  the 
jump  conditions  at  a  characteristic  surface 
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0 


(35) 


c  \\ip_]  +  i  \ii±]  .  *  \i*. il  +  / 

*/|**/J  *  \jx,'\  v  ]*»;]  i  L**,: 


lfx  it 

/>  *  dx, 


r.  *  (*** 

t  _ 


t  ir  fy  +  urf 


•I  [ft] 


7r-  M  *  f‘ 


u.-fx  +  ir/y  v-  ur-f^ 


l**/ 


fh  fe] + (»?*  *  ^ +  •**, )  [ff ]  ■ 0  (3,) 

*  u-fy  +  ur?f  j  [^4j  1  0  i,  ~ E+1,E+2,  -  ,NS  (40) 


v-  ia/  v-  or/’ 


'- )  [ft]  ■ 


«/  =  f,Z.-  -,fVV  (4!) 


The  vanishing  of  the  determinant  of  the  coefficients  (Equation  7)  implies  that 

/ 

for  X>t  -  0  to  be  a  characteristic  surface  the  condition 

( *  ^y  *  "fy  Y*  **"  3j(“'*  *  +***$-**(!*  ^  ^a)}  3  0  (42) 

must  be  satisfied.  It  can  be  seen  from  Equation  (42)that  there  are  two  sets 
of  characteristic  surfaces  corresponding  to  each  of  the  two  distinct  factors. 
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Each  factor  in (42)corresponds  to  conditions  for  the  directions 
of  the  normals  ,  fy  )  to  the  respective  characteristic  surfaces. 

The  first  factor  is  the  dot  product  of  the  velocity  vector 


>’ 


U>Z  +  V*j  *  Uf'j. 


(43) 


and  the  vector  normal  to  the  characteristic  surfaces  -P(X.  )  =  O 


vf  =  ~  r  * 

3p 


(44) 


Hence,  the  vanishing  of  this  product  implies  that  the  surfaces  composed  of 
streamlines  are  bJS-E  +  NV+  2  -  fold  characteristics  corresponding  to  the 
power  of  the  term.  This  fact  will  prove  to  be  very  useful  in  carrying  out 
a  numerical  integration  scheme  for  the  full  nonequilibrium  equations. 

The  second  factor  of (42)implies  thct  the  normal  vector 
(■fj,  -fy  ,  -f^)  to  the  surface 

=  constant  (45) 


lies  on  a  cone  of  the  second  degree  given  by 

a,ikxix'k  -  urx3')Z  ~  a,2(x’*4-  x/2)  =  0  (46) 

Equation  (46)  can  be  shown  to  be  real  and  thus  the  governing  equations  are 
hyperbolic  whenever  (U*+  ts  1  +  ur  1  )>&  i.  e.  when  the  flow  is  supersonic. 

Any  displacement  (  dx ,  dLyy  df.)  along  the  characteristic 
surface  (45)  through  a  point  P  can  be  given  in  the  limit  by 


f^dx  +  f^dy  +  ?yd$  =  0  (47) 

and  it  follows  that  if  ( -fx  ,  Ty  ,  )  lies  on  the  normal  cone  (46)  through 

P,  then  the  envelope  of  all  such  displacements  through  P  is  itself  a 
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quadratic  cone.  This  envelope,  which  ;s  called  the  "characteristic  conoid" 
at  P  is  given  by  (See  Courant  and  Hilbert  P.  563)^  inverting  the  quadratic 
form  (46)  as 

eb-j~' *  (tr*+  ur2 - a*)dj}M  +  (u>M+ur2-atl)dt/2+(ui-tir,-aa)ci^i 

(48) 

-  2 u,  t/'cL# dy  -  2uuroCxoL/  -  2 is* us* at y at j  -  O 

where  is  the  inverse  of  (X-j  in  (46).  Equation  (48)  can  also  be 

shown  to  be  the  equation  of  the  well  known  Mach  conoid  which  usually 
appears  in  supersonic  flow.  In  general,  it  represents  a  curvilinear  cone, 
whose  coefficients  are  functions  of  the  coordinates  (X,y,}  )  but  whose 
generators  are  tangent  to  the  generators  of  the  local  Mach  cone  at  F.  In 
any  numerical  procedure,  the  characteristic  conoid  is  replaced  by  an 
average  approximation  to  the  local  Mach  cone. 

The  generators  dX  •  (called  the  bicharacteristics)  of  the  local 

2  1 

Mach  cone  (48)  can  be  expressed  parametrically  by 

dXl  -  (a-  +  Cat i  cos 6  *  C/S-  sin 6)  dr  (49) 


where  S  is  a  parametric  angle,  ,  JS>^ 
C  is  defined  by 

C‘.  J!f- 

c  f '-a'  m‘ 


are  an  orthogonal  set  of  vectors, 


77  (50) 


and  7  is  the  time  taken  for  propagation  of  a  wave  along  a  ray  corresponding 
to  a  bicharacteristic. 

To  determine  the  cct  and  jSt  ,  we  now  consider  a  special 
intrinsic  coordinate  system  of  the  type  originally  proposed  by  Thornhill 
(Figure  1).  Here,  <f  is  the  angle  made  by  the  plane  containing  the  lines  j  and 
i  with  the  plane  containing  the  lines  y  and  y  .  L  is  a  ray  of  the  character¬ 
istic  conic  (Mach  surface)  making  an  angle  yd  with  the  velocity  vector  y  and 
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G,  p  are  the  polar  angles  of  the  velocity  vector.  The  coordinate  system 
is  chosen  aligned  as  shown  with  £  along  ^  ,  7j  in  the  pl,  ne  of  y  and 
y  perpendicular  to  ^  and  Q  perpendicular  to  the  plane  of  y  and  ^  . 
i*t  and  %  are  chosen  to  be  unit  vectors  and  are  given  by 


cos  0  cos  pT  +  sin  &J  +  cos  6  sin  p  K  *■  ~  (ul+  wj  +  ur~k  ) 


rj  =  -Sin  9  cos  pi  -hcosdj  -Stn&sm  pk  -  afI+  +  CZ3  k  (52) 


^  *  -sin  pi  +  cos  pi  *  i  *  j  -t-  E  (53) 

Therefore,  in  the  form  of  (49)  the  bicharacteristic  direction  Z  is  given 
by 

dx  -  (u,  -  C  sin  9  cos  p  cos  &  -  C  sm  p  sinS)dT 

dy  =  (tr+  C  cos  9  cos &)dr  (54) 

dj.  =  (ur  -  Cstn  p) suo  0cosS  +  Ceos  p stn  S )dr 

2.4  THE  COMPATIBILITY  EQUATIONS 

As  shown  in  Paragraph  2.  1,  the  compatibility  equations  are 
determined  by  forming  a  linear  combination  of  the  flow  equations  and 
choosing  the  coefficients  such  that  the  derivatives  in  the  direction  normal 
of  the  characteristic  surface  disappear.  There  are  two  sets  of  compat¬ 
ibility  equations  corresponding  to  the  two  sets  of  characteristic  surfaces 
which  were  described  in  the  last  section.  Actually  since  the  streamlines 
are  e/v'Pi -fold  characteristics,  there  are  NS  -E+N/+-  %  compatibility 

equations  which  can  be  found  to  apply  along  them. 

The  equations  which  apply  along  the  streamline  surfaces  could 
be  found  by  applying  the  method  developed  in  2.1  to  the  full  equations  of 
motion.  However,  in  this  case,  it  is  easier  to  recall  that  the  equations 
sought  do  not  contain  derivatives  in  the  direction  normal  to  the  streamline. 
This  property  can  be  used  to  choose  the  required  relations  by  inspection 
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0 


from  the  original  partial  differential  equations.  The  streamline  directions 
are  give  n  by 


di  -  $dt  (55) 

where  §  is  the  coordinate  along  the  streamline.  Thus,  MV-£  *  A /S 
compatibility  equations  can  be  found  by  transforming  Equations  (17)  and  (19) 
to  directional  derivatives  along  the  streamline.  The  equations  become 


and 


im  Eft,  •  •  •,  NS  ( 56) 


J-  t,2.-  ■  ‘.NV  (57) 


By  introducing  Equation  (20)  for  the  enthalpy/?  into  the  energy  equation  (16) 
and  differentiating  in  the  manner  used  to  obtain  Equation  (29),  the  energy 
quotation  along  a  streamline  can  be  written  as 


A.  dp 
7  2?* 


=■  o 


(58) 


One  further  relation  is  required.  This  is  found  from  scalar  multiplication 
of  the  momentum  equations  and  the  velocity  vector  5  as 


d  /_/  i  j  _  _  J_  /dp  \ 

df  \2  ?  J  ~  \d4  ) 


£3 


J_  dh 
A 


(59) 


where  the  last  term  on  the  right  hand  side  of  (59)  is  obtained  from  the 
energy  equation  (16).  Since  Equations  (56)  -  (59)  contain  no  derivatives 
in  the  direction  normal  to  the  streamlines,  they  are  the  required 
Nd-£  +NV+2 compatibility  equations. 
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The  compatibility  equation  which  applies  along  the  characteristic 
conoid  (48)  can  be  found  by  applying  the  method  of  Paragraph  2.  1  to 

Equations  (13)  and  (32).  Equations  (35)  -  (39)  are  representative  of  the 

/  /  • 

equations  of  motion  transformed  to  characteristic  coordinates  (  )• 

Multiplying  these  equations  by  -^j^^and  respectively  leads  to  the  following 
equation  for  the  term  flM  p  (f)  presented  in  Paragraph  2.  1. 


^  (-p *,  (f  vfhuj. 

*  (h/i-Tt)  *  h,  (y) *(V (I  V?)  *  M J  [33^  *  W  <*  ^ 


(60) 


The  condition  for  which  the  surface  })  -  0  is  to  be  a  characteristic 

surface  is  that 


V  )  *  O  (61) 

Hence,  the  coefficients  of  must  all  vanish  identically.  Equation  (60) 

thus  leads  to  the  following  expressions  for  the  linear  factors  h „  . 


Using  hese  values  for  h ,  the  compatibility  equation  corresponding  to 
Equation  (11)  become s 


■ 

where  <p  ~^,h  for  l  =  2,3  respectively. 
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Equation  (65)  :an  be  simplified  by  choosing,  as  the  coordinates 
(  X/1'J  )  the  basic  characteristic  coordinate  system.  Referring 

again  to  Figure  1,  take  as  the  coordinates  ,  /V  and  87  where  L 
is  along  the  bicharacteristic  given  by  S  ,  V  is  orthogonal  to  L  in  the 
plane  of  fj  and  £  and  id  =  Lx  /7  is  orthogonal  to  the  characteristic 
surface  and  positive  outward.  Here,(/^',  X>^' 3  &9')  =*  (  M ,L,  /V) 
and 

■fi  *  M  s  cos p  cos &rj  +  cos  J5  sin  SC  ~  Sen  /SC  (66) 

»  </'  *  L  -  sen p  cos  S7j  ■+•  sin p sen  6  £  +  cos  p/[  (67) 


h  11  N  »  sin  Srj  -  cos  £  C 


where  A  =  Sen  — — 
'  A/ 


With  respect  to  the  characteristic  coordinate  system  (46)  -  (48)  Equation 
(65)  becomes 

/  & p  /  2 

—  — -  +  (- sinp  cosS  cos  0 cos  /  -  sin/Osin 6 sen / -cos  p cos S sen 6 cos d> /sen £ 

pa>  oL  r/ 

-cos/SsmSsin //s*np)jjr  +  (cos 6  sen  </ -sin  8 Sen 0cos  </) 

4- (sin P cos S cos  8  v»  cos/? cos  8  cos  0/sinA)^X  ■+•  sen  6 cos  0 

(70) 

+  (sinfl  sen  £cos  -  cos  6  sen  &  Sen  </  senfi 

-  cos*p  cos  £  send  sen  <//stn/3  +  cos  Sen  S  cos  /e/sinp 

d L 

+  ( -sin  8 sin  8 sin  </>  -  cosS cos  </)  ~  * 

Changing  from  the  variables  (  p,  u,  ir,  un )  to  the  variables  ( p,  ) 

and  simplifying,  the  final  form  of  the  compatibility  equation  becomes 

co  ip  dp  „  38  „ 

- j-  — r  *  cosS-—-  +  cos  0  sen£  — — 

/sg  01  dL  dl  /71) 


where  _£?  and  — — 


/  30  \ 

=  -  sen p  [cos  S  cos  0 —  -senf>-^  -f  ?  J 

are  the  derivatives  along  and  normal  to  the 
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bicharacteristic  6  .  Note  that  Equation  (71)  contains  derivatives  of 

P‘  &  an<3  only.  Equations  (56),  (57),  (58)  and  (71)  are  the  basic  set  of 
equations  which  will  ,  ised  to  develop  a  numerical  integration  scheme  for 
three-dimensional  steady  flow  in  Section  3. 
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Section  3 

THE  NUMERICAL  SOLUTION 

In  this  section,  the  details  of  a  multivariable  numerical  integration 
scheme,  which  is  based  upon  the  method  of  characteristics  applied  to  the 
equations  developed  in  Section  2,  are  presented.  First,  the  general  require¬ 
ments  for  the  stability  and  convergence  of  such  a  scheme  are  considered, 
and  then  several  possible  difference  networks  are  discussed  and  compared. 
Finally,  the  details  of  the  proposed  techniques  are  given  for  both  steady  and 
nonsteady  nonequilibrium  flows. 

3.  1  STABILITY  AND  COVERGENCF 

In  ail  numerical  approaches,  the  solutions  to  a  system  of  partial 
differential  equations  are  represented  by  values  of  the  dependent  variables 
for  certain  discrete  values  of  the  independent  variables.  In  general,  using 
this  procedure,  there  are  three  types  of  errors  which  would  cause  the 
numerical  values  to  be  different  from  those  which  would  be  obtained  from 
the  exact  solution  of  the  partial  differential  equations  themselves.  These 
can  be  classified  as  truncation  errors,  round  off  errors,  and  errors  which 
appear  in  the  initial  values.  The  first  is  the  error  which  results  from 
using  a  difference  formula  which  is  an  approximation  to  the  true  equation 
and  may  be  considered  as  the  error  incurred  from  representing  an  infinite 
series  by  a  finite  number  of  terms.  The  second  arises  from  the  need  to 
use  finite  decimal  numbers  in  the  computation,  while  the  last  may  occur 
because  of  inaccurate  description  of  the  initial  and  boundary  conditions. 

In  principle,  the  cumulative  effects  of  these  errors  upon  the  solution  can 
be  kept  under  control  by  reducing  the  mesh  spacing  between  grid  points  and 
carrying  out  the  calculations  with  sufficiently  high  precision. 

In  practice,  however,  the  cumulative  departure  from  the  exact 

solu.  on  for  a  fixed  range  of  the  independent  variable  will  usually  grow 

unboundedly  as  the  mesh  size  approaches  zero.  This  is  because  the  precision 

required  to  sufficiently  control  the  growth  would  far  exceed  the  capacity  of 

existing  computers.  If  this  growth  occurs  exponentially  as  the  mesh  size 

decreases,  it  is  generally  considered  unmanageable  for  computational 

7 

purposes,  and  the  procedure  is  termed  unstable.  The  establishment  of 
conditions  necessary  to  assure  stable  solutions  thus  assumes  a  major  role 
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in  any  numerical  scheme.  For  a  set  of  first-order  linear  partial  differential 

40 

equations  in  n-independent  variables,  Courant,  Friedrichs,  and  Lewy,  in 

their  classical  paper,  have  shown  that  a  necessary  conditions  for  stability 

(i.  e.  ,  the  C-F-L  stability  condition)  is  that  the  domain  of  dependence  of  the 

partial  differential  equations  (e.  g.  ,  the  circle  £  + 'J*  ~  for  the 

Mach  cone  in  steady  flow  shown  in  Figure  1)  be  contained  within  the  domain 

of  dependence  of  the  difference  equations.  This  latter  domain  is  usually 

called  the  "convex  hull"  of  the  difference  scheme  and  is  the  region  obtained 

by  connecting  with  straight-line  segments  the  outermost  points  used  in 

7 

solving  the  difference  equations.  Forsythe  indicates  that,  for  a  central 

differencing  scheme  using  four  initial  points  in  three  dimensions,  or  a 

comparable  number  in  n-di.mensi  ons,  the  C-F-L  conditions  is  both  necessary 

and  sufficient  to  guarantee  stability,  and  thus  also  convergence  of  the  differ- 
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ence  scheme.  Hohn  has  also  found  this  to  be  the  case  when  simplical 
difference  schemes  (i.  e.  ,  schemes  where  a  minimum  number  of  points  are 
utilized  in  the  initial  plane)  are  used  for  solving  linear  equations. 

While  the  equations  of  motion  presented  in  Section  2.  2  are  not 
linear,  it  can  be  argued  on  physical  grounds  that  the  C-F-L  stability  criterion 
is  also  likely  necessary  for  stability  of  a  difference  scheme  based  upon 
these  equations.  Consider,  for  instance,  the  simple  case  of  supersonic 
flow  parallel  to  the  x-axis  and  assume  that  the  four  points  (  0,  +  h  ,  0), 

(0,-  hf  ,  0)  ,  ( 0,  0,  -  ht  )  ,  and  (0,0,  hz)  comprise  the  points  in  the 
initial  plane  X  =  0  to  be  used  in  a  central  differencing  scheme  for  calcu¬ 
lating  the  flow  properties  at  a  new  point  (k.,  0,0),  (See  Figure  2).  The 
domain  of  dependence  of  the  difference  equations  is  then  the  rhombus 
formed  with  the  four  points  in  the  initial  plane  as  vertices.  If  the  C-F-L 
stability  condition  is  not  satisfied,  the  circle  which  is  the  intersection  of 
the  Mach  forecone  f-om  (k.,0,0)  with  the  initial  plane  (e.  g.  ,  the  domain 
of  dependence  of  the  differential  equations)  would  not  be  contained  within 
the.  rhombus.  Hence,  certain  initial  values  (e.  g.  ,  points  within  the  circle 
but  outside  the  rhombus)  could  be  arbitrarily  changed  without  affecting  the 
solution  of  the  difference  equations.  The  departure  of  the  difference  solution 
from  the  solution  of  the  differential  equations  could  thus  become  very  large 
and  the  scheme  would  be  unstable. 
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Mathematically,  Strang  '  has  shown  that  the  convergence  of  the 
nonlinear  equations  depends  uniquely  upon  the  stability  of  the  linear  equations 
for  cases  where  the  equations  and  their  solution  possess  enough  continuous 
derivatives.  For  many  cases  of  fluid  flow,  where  the  gradients  are  not  too 
high  (e.  g.  ,  excluding  flow  throueh  shocks,  etc.)  these  conditions  would  be 
expected  to  be  met,  and  the  C-F-L  condition  should  be  sufficient  to  insure 
stability  using  the  simple  differencing  schemes  discussed  above. 

Although  the  above  discussion  applies  specifically  to  numerical 
solutions  using  a  direct  finite  difference  network,  the  physical  arguments 
given  in  paragraph  three  can  also  be  seen  to  hold  when  a  method  of  char¬ 
acteristics  approach  is  used.  Now,  the  four  vertex  points  of  the  rhombus 
(convex  hull  of  difference  equations)  would  represent  base  points  of  four 
bicharacteristics  along  which  the  equations  are  to  be  integrated.  These 
po;nts  would  now  lie  on  the  circle  and  the  rhombus  would  thus  be  contained 
completely  within  the  circle.  Hence,  when  no  other  points  are  included, 
the  scheme  would  necessarily  be  unstable  because  some  of  the  initial  data 
is  neglected.  In  the  technique  described  below,  the  over-all  procedure 
utilizing  a  method-of-characteristics  solution  is  stabilized  by  introducing 
initial  data  points  outside  the  domain  of  dependence  of  the  differential  equa¬ 
tions  and  obtaining  the  flow  properties  at  the  bicharacteristic  base  points 

by  interpolation  among  these  initial  points  along  the  lines  employed  by 
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Sauerwein  and  Sussman,  Tsung,  “  and  Butler  and  Talbot.  The  extra 

points  added  now  form  the  outer  boundary  of  the  convex  hull  of  the  charac¬ 
teristic  difference  network  and  the  C  -F-L  stability  condition  can  be  modified 
to  apply  to  this  region. 

3.  2  CHOICE  OF  A  PRATICAL  INTEGRATION  NETWORK 

The  method  of  characteristics  has  been  applied  extensively  to 
solve  two-dimensional  flow-field  problems.  In  this  case,  there  are  two 
characteristics  which  pass  through  any  initial  point.  The  compatibility 
equations  are  ordinary  differential  equations  which  can  be  written  in 
finite  difference  form  and  used  to  solve  for  the  flow  variables  at  a  new 
point.  Two  basic  networks  have  been  utilized  for  carrying  out  the  numerical 

1  z 

integration.  In  the  first,  as  explained  by  Ferri  u  and  shown  in  Figure  3 
point  is  located  as  the  intersection  of  two  characteristics  from  initial 
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Figure  2  COURANT-FRI EDRICHS-LEWY  CONDITION 
FOR  FINITE  DIFFERENCE  SCHEME 


Figure  3  FORWARD  CHARACTERISTIC  NET  FOR 
TWO  INDEPENDENT  VARIABLES 
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points  Pj  and  P  ,  respectively.  The  flow  variables  at  point  P^  are  obtained 
by  simultaneously  solving  the  compatibility  equations  written  along  these 
characteristics.  Here,  compatibility  equations  along  the  streamline  P^  P^ 
may  be  introduced  as  auxiliary  conditions.  The  second  method  shown  in 
Figure  4  was  originally  proposed  by  Hartree.  This  consists  first  of  sel¬ 
ecting  the  point  P^  off  the  initial  line  and  intersecting  the  characteristics 
through  P  with  the  initial  line  at  points  P^  and  P^.  The  flow  properties 
at  points  P^  and  P7  are  determined  by  interpolation  and  the  finite  difference 
forms  of  the  bicharacteristic  equation  are  used  to  solve  for  the  flow  vari¬ 
ables  at  point  P^.  Here  again,  the  streamline  through  P^  can  be  supplied 
when  more  than  two  dependent  variables  are  involved. 

Both  methods  describee  above  for  the  two  variable  problem  are 
essentially  the  same,  because  regardless  of  how  point  P^  is  located,  the 
choice  of  characteristic  directions  through  P^  is  not  arbitrary.  However, 
this  is  not  the  case  for  problems  involving  more  than  two  independent 
variables.  Indeed,  as  was  shown  in  Equation  (49),  the  generators 
of  the  characteristic  cones  in  three -independent  variables  comprise  a 
one -parameter  family  of  curves  through  any  point  P.  Thus,  there  is  a 
certain  degree  of  arbitrariness  involved  in  choosing  a  basic  finite  differ¬ 
ence  network  for  use  in  a  multivariable  characteristic  procedure. 

Before  discussing  some  of  the  networks  which  have  thus  far 
been  proposed,  it  will  first  be  useful  to  introduce  the  generalization  of  the 
two  characteristic  directions  associated  with  any  point  P  in  two  dimensions. 
For  this  purpose,  consider  an  arbitrary  curve  that  acts  as  a  source  of 
disturbance  in  three-space,  (see  Figure  5).  A  characteristic  conoid  is 
associated  with  each  point  on  the  curve  in  the  manner  indicated.  The 
envelope  of  all  these  conoids  .which  also  contains  the  original  curve,  forms 
two  distinct  surfaces  (  and  Z.2  )  which,  since  they  are  comprised  of 
bicharacteristics,  are  also  characteristic.  Hence,  through  any  curve  in 
three  space,  there  are  two  characteristic  surfaces  corresponding  to  the 
two  characteristic  lines  through  a  point  in  two  space.  The  geometry 
involved  in  the  three-dimensional  method  of  characteristics  has  thus  become 
evident.  Through  each  point  P  in  space,  the  normal  N  to  a  characteristic 
surfact  lies  on  the  normal  cone  given  by  Equation  (45)  (see  Figure  6). 
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Figure  4  BACKWARD  CHARACTERISTIC  NET 
TWO  INDEPENDENT  VARIABLES 


Figure  5  CHARACTERISTIC  SURFACES  THROUGH  ARBITRARY  CURVE 
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The  envelope  of  all  such  surfaces  through  P  is  just  the  characteristic  conoid 
at  P  and  is  itself  a  characteristic  surface.  The  curves  of  contact  between 
the  characteristic  conoid  and  the  enveloping  surfaces  are  the  bicharacteristic 
curves.  Thus,  the  characteristic  conoid  may  be  regarded  as  being  generated 
by  the  bicharacteristic  curves  through  P,  its  vertex.  This  characteristic 
conoid  coincides  with  the  Mach  conoid  at  every  point  on  its  surface.  As  shown 
in  Figure  7,  the  boundary  of  the  Mach  conoid  is  just  the  envelope  of  all  the 
local  infinitesimal  Mach  cones.  The  characteristic  conoid  is  thus  everywhere 
tangent  to  a  local  Mach  cone. 

The  solution  of  flow  field  problems  by  the  method  of  character¬ 
istics  in  more  than  two  independent  variables  consists,  therefore,  of 
choosing  from  the  infinity  of  directions  available  a  particular  set  of 
bicharacteristics  or  characteristic  surfaces,  writing  the  equations  devel¬ 
oped  in  Section  2  along  them  in  finite  difference  form  and  solving  the 
resulting  numerical  equations  simultaneously.  Various  finite  difference 

networks  have  been  proposed  for  accomplishing  these  integrations.  The 
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most  important  of  the.  ;  have  been  discussed  and  compared  previously.  ’ 

Each  of  the  proposed  schemes  are  summarized  here  for  ease  of  reference 
and  finally  evaluated  from  the  standpoint  of  their  applicability  to  finite-rate 
chemically  reacting  flow  field  studies.  In  ali  of  the  procedures,  the  char¬ 
acteristic  surfaces  and  conoids  discussed  above  are  replaced  by  suitable 
average  planes  and  cones.  The  nomenclature  used  by  Fowell  and  Sauerwein 
is  repeated  here. 

1 4 

Tho  rnhill  originally  proposed  two  integration  schemes.  In  the 

first,  termed  the  "tetrahedral  bicharacteristic  line  network"  (Figure  8),  the 
new  point  P^  is  located  as  the  common  intersection  point  of  Mach  cones  from 
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Figure  6  NORMAL  CONE  AND  CHARACTERISTIC  CONE 


Figure  7  MACH  CONOID  AS  ENVELOPE  OF  LOCAL  MACH  CONES 
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each  of  three  initial  points  P^,  P^,  and  P^.  The  li:  es  P^  P^,  P^  P^  and 
Pj  P^  then  represent  numerical  approximations  to  three  bicharacteristics 
and  the  streamline  from  point  P^  P^  to  the  initial  plane  provides  a  fourth 
bicharacteristic  direction.  The  properties  at  Pg  are  determined  by  inter¬ 
polation.  The  equations  developed  in  Section  2  written  in  finite  difference 
form  along  these  bicharacteristic  directions  are  used  to  solve  for  the  flow 
properties  at  the  Point  P  .  The  main  advantage  of  this  network  was  the  fact 
that  the  base  points  P  ,  P?)  and  P  remained  fixed  throughout  the  iteration 
process.  However,  Sauerwein,  7  in  attempting  to  apply  the  technique  to  non¬ 
steady  flow  problems,  determined  that  it  was  numerically  unstable.  The 
reason  for  this,  as  shown  in  Figure  8,  is  that  the  domain  of  dependence  of  the 

differential  equations  (the  circle  P^  P^  P^)  is  not  contained  within  the  domain 
of  dependence  of  the  difference  equations  (the  triangle  P^  P,)  and  the 
C-F-L  stability  condition  as  introduced  in  Section  3.  1  is  continually  vio¬ 
lated.  Sauerwein  then  proposed  what  he  termed  the  "modified  tetrahedral 
characteristic  line  network.  "  In  this  network  (Figure  9),  three  points 
F* i ^ »  ^23’  anc*  representing  the  points  of  tangency  of  the  circle  in¬ 
scribed  within  the  original  triangle  P^P^P^  are  chosen  as  new  initial 
points.  The  properties  at  points  P^,  ^*23’  anc*  ^31  are  determined  by 
interpolation  and  the  method  then  follow's  the  procedure  described  above 
for  the  unmodified  network.  The  C-F-L  condition  for  the  points  P^,  ^*23’ 
and  P^  is  seen  to  be  satisfied  and  the  technique  was  found  to  be  stable. 

The  second  method  proposed  by  Thornhill  is  the  "tetrahedral 
bicharacteristic  surface  network"  (Figure  10).  Again,  three  noncolinear 
points  Pj,  P^i  and  P^  are  chosen  in  the  initial  plane.  Point  P^  is  deter¬ 
mined  as  the  mutual  intersection  point  of  the  three  average  internal 
characteristic  planes  through  lines  P  P?,  P  P  and  P  P  .  The  Mach 
forecone  from  P.  defines  three  bicharacteristic  lines  P  .  Pr,  P,  P. ,  P,  P„ 
as  the  lines  of  tangency  between  the  cone  and  each  of  the  original  charac¬ 
teristic  planes.  The  streamline  P^  Pg  again  furnishes  the  fourth  bichar¬ 
acteristic  direction.  The  flow  properties  at  the  base  points  Pc,  P , ,  P_ 
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and  Pg  are  determined  by  interpolation  and  the  flow  properties  of  point  P^ 
are  then  solved  for  along  the  bicharacteristics  as  before.  This  method  and 
the  modified  method  discussed  above  are  essentially  the  same,  as  the 
bicharacteristics  chosen  run  from  the  tangent  points  of  the  circle  inscribed 
in  an  initial  triangle  to  the  new  point  at  which  the  flow  properties  are  desired. 
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‘-CIRCLE  OF  DEPENDENCE 
OF  DIFFERENTIAL  EQUATIONS 


Figure  8  TETRAHEDRAL  B I  CHARACTER  I  ST  1C  LINE  NETWORK 


Figure  9  MODIFIED  TETRAHEDRAL  BICHARACTERISTIC  LINE  NETWORK 

36 


Hence,  since  the  "modified"  characteristic  network  led  to  a  stable  solution, 
it  seems  plausible  to  suspect  that  the  present  scheme  would  also  be  stable. 
Indeed,  such  is  the  case,  as  was  shown  by  Tsung  in  applying  this  technique 
to  the  simplified  case  of  three-dimensional  potential  flow  problems. 

A  "network  of  intersections  of  reference  planes  with  character- 

18  43 

istic  surfaces"  (Figure  11)  as  discussed  first  by  Ferrari,  Sauer,  and 
44 

Ferri  considers  two  sets  of  orthogonal  coordinate  planes  (e.  g.  ,  y  =  constant 
and /C  =  constant)  and  utilizes  a  two-dimensional  characteristics  network. 

The  intersection  of  the  two  coordinate  planes  define  lines  and  Pg  P^ 

from  which  the  integration  proceeds.  Points  Pg  and  Pg  are  located  by  a 
forward  marching  two-dimensional  characteristics  scheme  whereby  the 
characteristics  are  confined  to  a  =  constant  planes.  Projections  of  the 
streamlines  through  points  Pg  and  Pg  on  to  the  x.  -  constant  planes  locate 
points  P^  and  Pg.  Lines  P^  Pg  and  Pg  Pg  then  serve  as  two-dimensional 
streamlines  which  can  be  used  as  auxiliary  conditions  for  determining  the 
flow  properties  at  points  Pg  and  Pg.  The  variations  of  properties  in  direc¬ 
tions  across  the  reference  planes  x  =  constant  are  obtained  by  using  poly¬ 
nomial  fits  to  evaluate  the  c ross -derivative  terms  which  appear  in  the 
equations  written  along  the  lines  Pj  Pg  and  P^  Pg  or  the  lines  Pg  Pg  and 

and  P.  P, .  In  order  to  evaluate  the  cross-derivatives  which  occur  in  the 
4  6 

equations,  data  points  must  be  given  in  planes  y  =  constant  for  each  step 

of  the  calculation.  Generally,  however,  the  newly  calculated  points  Pg  and 

Pg  wil1  not  lie  in  the  same  }  -  constant  plane.  Interpolation  or  extrapolation 

to  point  Pg  in  the  desired  ^  =  constant  plane  is  thus  required.  While 

the  integration  network  leading  to  the  new  field  points  in  this  method  does 
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not  consist  of  bicharacteristic  curves,  Ferrari  has  shown  that  the 

equations  used  are  completely  equivalent  to  the  original  differential  equations 

and  that  the  existence  and  uniqueness  tests  devised  by  Titt  are  also  satisfied 

2  3 

here.  Morretti,  et  al  have,  in  fact,  obtained  practical  results  using  this 
scheme.  li  v'ould  appear  that  the  reason  Morretti  was  able  to  obtain  stable 
results  is-  that  the  points  used  as  base  points  for  the  two-dimensional  charac¬ 
teristics  plus  the  points  used  to  determine  the  c ross -derivatives  all  lie  out¬ 
side  the  Mach  forecone  from  the  new  point.  Whenever  flow  conditions  are 
such  that  the  points  used  do  not  include  the  Mach  forecone,  the  procedure 
would  probably  be  unstable  because  the  C-F-L  conditions  would  be  continually 
violated.  Another  disadvantage  of  this  network  is  that  interpolation  must  be 
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performed  throughout  the  flow  field  in  order  to  confine  the  newly  calculated 
points  to  planes  of  ^  =  constant.  Usually,  linear  interpolations  are  used 

and  the  over-all  effects  of  these  linear  interpolations  can  appreciably  affect 
the  accuracy  of  the  technique. 

A  fifth  technique  which  has  had  some  practical  application  is 
21 

that  proposed  by  Butler,  (Figure  12).  This  procedure  represents  an 
extension  of  Hartree's  method  as  described  above  for  the  two-independent 
variable  problem.  Given  the  initial  data  plane,  a  new  point  P,.  is  chosen 
in  a  new  data  plane.  The  coordinates  of  are  thus  known  exactly.  A 
first  estimation  of  the  flow  variables  at  P,.  is  made  and  four  bicharacteristics 
(or  one  more  than  the  minimum  required)  located  in  planes  90®  to  each  other 
are  passed  back  to  the  initial  surface.  Linear  combinations  of  the  equations 
along  these  four  bicharacteristics  result  in  relations  involving  derivatives 
in  the  bicharacteristic  directions  only  at  the  new  point.  The  streamline 
P^  P,.  from  Pj-  provides  another  relation  which  is  required  for  obtaining 
a  complete  solution.  T..e  flow  properties  at  the  base  points  Pj,  P^,  P^»  P^, 
and  are  determined  by  interpolation.  Butler  and  Talbot  have  succeeded 
in  applying  this  technique  to  some  simple  two-dimensional  unsteady  flow 
problems.  A  searching  scheme  was  used  to  determine  the  nine  points  which 
were  closest  to  the  particular  base  point,  say  Pj,  being  fit.  Then  bivariate 
Lagrange  interpolation  through  these  nine  points  was  used  to  determine  the 
flow  properties  at  P^.  However,  the  procedure  was  found  to  be  stable  only 
when  all  of  the  base  points  associated  with  a  given  calculation  were  confined 
within  the  same  nine-point  fit,  and  hence  only  when  the  C -F -L'condition  was 
satisfied.  If  one  set  of  points  is  used  to  fit  a  given  base  point,  while  another 
set  of  points  is  used  to  fit  the  remaining  base  points,  the  procedure  was  found 
to  be  unstable.  This  can  probably  be  traced  to  the  fact  that  not  all  the  points 
in  the  domain  of  dependence  of  the  difference  equation  are  treated  with  equal 
weight  and  hence  the  C-F-L  condition  is  in  principle  violated 

Other  methods  have  been  proposed  for  which  no  practical  three- 
dimensional  results  have  as  yet  been  obtained.  Holt  using  the  work  of 
Coburn  and  Dolph,  introduced  a  network  which  Fcwell  called  a  prismatic 
network  of  characteristic  surfaces,  (see  Figure  13).  This  method,  which 
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Figure  12  PENTAHEDRAL  BICHARACTERISTIC  LINE  NETWORK 
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is  similar  to  the  network  of  intersection  of  reference  planes  with 
characteristic  surfaces  was  devised  because  of  the  misgivings  Coburn 
and  Dolph  had  regarding  the  use  of  non -bicharacte ristic  directions  along 
which  the  integrations  were  to  be  performed.  As  mentioned  above, 
however,  Ferrari  later  showed  that  this  objection  was  clearly  unjustified. 
The  method  consists  of  choosing  the  two  bicharacteristic  directions  P^P,. 
and  P,P  as  two  coordinate  directions  at  a  general  point.  The  third 
direction  is  provided  by  the  line  P^P^  initial  plane.  In  general,  this 

method  requires  that  the  flow  properties  be  known  on  some  surface  P^P^P, 
other  than  the  initial  surface  (e.  g.  ,  plane  of  symmetry,  etc.)  and  this 
would  seem  to  limit  its  application. 

These  are  essentially  all  of  the  basic  characteristic  nets  which 

have  thus  far  been  advanced.  Slight  variations  of  the  method  of  intersection 

of  reference  planes  and  characteristic  surfaces  have  been  discussed  by 
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Sauer  and  Shaetz  (called  by  them  the  method  of  "near  characteristics") 

24 

for  unsteady  flow  and  by  Kackova  and  Cuskin  for  steady  ilow.  These 
methods  would  all  seem  to  have  the  same  advantages  and  disadvantages  oi 
the  Ferrari,  Sauer,  Ferri  technique  descmbed  above. 

All  of  the  above  integration  procedures  are  seen  to  require  some 

sort  of  interpolation  at  the  base  of  the  bicharacteristics  in  the  initial  data 

surface.  The  "tetrahedral  bicharacteristic  line  network"  which  sought  to 

avoid  this  problem  was  found  to  be  unstable.  Also,  for  every  method, 

interpolations  are  required  at  the  base  of  the  streamline  in  the  initial  plane. 

A  critical  point  to  be  considered  in  a  practical  application  of  any  of  the 

above  techniques  then  is  how  to  keep  the  number  of  interpolations  to  a 

minimum.  This  is  necessary  both  from  the  point  of  view  of  conserving 

machine  time  and  obtaining  accurate  representations  of  flow  fields.  When 

nonequilibrium  chemistry  is  included,  this  point  becomes  even  mere  salient. 

Here  a  total  number  of  A/S~£  interpolations  would  be  necessary  for 

the  compositions  and  energies  at  the  base  of  each  streamline.  However, 
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Curtis  has  shown  that  interpolations  for  species  which  vary  rapidly  across 
the  shock  layer  are  a  constant  source  of  trouble  leading  to  oscillations  in 
species  concc. Aiations  along  streamlines.  If  possible,  such  interpolations 
should  therefore  be  avoided. 
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Other  considerations  should  also  be  taken  into  account  when 

finite -rate  calculations  are  desired.  First,  experience  has  shown  that 

the  concentration  gradients  are  usually  much  steeper  along  streamlines 

than  are  the  gradients  in  other  flow  quantities.  Close  control  of  integration 

step  size  along  streamlines  is  thus  required  in  order  to  obtain  reasonably 

28  34 

correct  concentrations.  Second,  Sedney  and  Wood  et.  al.  have  shown 
that  the  use  of  entropy  as  a  dependent  variable  leads  to  large  errors  and 
should  be  avoided.  The  energy  equation  as  developed  in  Paragraph  2.  2  has 
avoided  this  problem  by  using  temperature  as  the  dependent  variable. 

Third,  since  any  finite-rate  chemistry  calculations  require  considerable 
machine  run  time  per  point,  a  method  of  reducing  the  number  of  points 
required  to  obtain  an  adequate  flow  field  representation  is  desired.  Probably 
the  best  way  to  reduce  the  number  of  calculations  is  to  keep  the  mesh  spacing 
as  uniform  as  possible. 

When  the  integration  techniques  discussed  above  are  assessed 
as  to  their  applicability  to  finite-rate  chemistry  calculations,  the  outlook  for 
any  of  the  networks  is  not  very  promising.  All  of  the  schemes  involve  inter¬ 
polations  for  concentrations  in  the  initial  plane.  The  network  of  intersections 
of  reference  planes  with  characteristic  surfaces  would  also  involve  inter¬ 
polations  on  the  concentrations  throughout  the  flow  field.  In  fact,  in 
order  to  control  mesh  size  spacings  by  either  adding  or  dropping  points 
all  of  the  methods  would  again  require  interpolations  throughout  the  entire 
flow  field.  The  methods  used  to  obtain  the  coordinates  of  the  new  point 
are  generally  quite  corpplex  in  all  cases  except  the  pentahedral  bicharacter- 
iatic  line  network  where  the  location  of  the  new  point  is  initially  given. 

Since  the  run  times  involved  in  finite-rate  chemistry  calculations  are  already 
long,  such  complicated  procedures  shf  ild  in  general  be  avoided. 

The  integration  network  described  here  seeks  to  avoid  the 

complicated  geometry  and  interpolation  problems  discussed  above.  It 

represents  an  extension  and  modification  of  the  two-variable  scheme 
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proposed  by  Curtis  and  essentially  combines  the  simplicity  of  the 
pentahedral  network  with  the  advantages  of  a  direct  marching  technique. 
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It  may  be  termed  a  "network  of  intersections  of  streamlines  with  reference 
planes"  since  the  basic  integration  coordinate  is  along  the  streamline  (or 
particle  lines  in  unsteady  flow)  and  new  points  are  located  on  successive 
reference  planes  or  surfaces. 

The  numerical  procedure  is  illustrated  for  a  single  field  point 

calculation  in  Figure  14.  All  the  flow  quantities - composition, 

pressure,  flow  angle,  velocity,  etc.  are  available  at  initial  data  points 
surrounding  point  Pj ,  in  the  non  characteristic  initial  plane  %  -  .  A  new 

plane  at  which  the  data  is  required  is  fixed  at  z  =  X^+  AX.  The  new  point 
is  first  located  as  the  intersection  of  the  streamline  from  Pj  with  the 
plane  x  =X^+AX,  Four  bicharacteristics  positioned  90“  apart  on  the  Mach 
forecone  from  P^  (Equation  54)  are  then  passed  back  to  the  initial 
plane  %  -  ,  The  flow  quantities  at  the  base  points  P^,  P^,  P,.  and  P^ 

are  determined  from  bivariate  surface  fitting  techniques,  remembering 
that  for  stability  the  region  fit  must  include  the  domain  of  dependence  of 
the  differential  equations.  Actually,  in  practice  the  surface  fit  for  each 
flow  variable  is  determined  only  once  during  the  iteration  process.  Thus, 
the  interpolations  become  merely  a  matter  of  evaluating  known  polynomial 
expressions.  The  velocity  ^  ,  temperature  T  .  concentrations  and 
vibrational  energies  6j  at  P^  are  found  by  integrating  Equations  (56) 
through  (59)  numerically  along  the  streamline  PjP^.  Pour  compatibility 
equations  (71)  written  in  finite  difference  form  along  the  four 
bicharacteristics  P^P^,  ^4^2’  ^*5^2  anc^  ^*6^2  appropriate  average 

values  of  the  remaining  flow  quantities  ---  pressure  p  ,  and  flow  angles 
6  and  .  This  procedure  is  then  iterated  to  convergence  of  the  flow 
properties  at  point  P^. 
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•  INITIAL  DATA  POINTS  FIT 
X  B I  CHARACTER  I  ST  I C  BASE  POINTS 
X  STREAMLINE  BASE  POINT 
O  NEW  DATA  POINTS 


Figure  14  PROPOSED  FIELD  POINT  NETWORK 
THREE  INDEPENDENT  VARIABLES 
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Perhaps  the  only  disadvantage  of  "the  network  of  intersections  of 
streamlines  and  reference  planes"  is  that  the  base  points  of  die  bicharac¬ 
teristics  tend  to  drift  around  in  the  initial  plane  during  the  iteration  cycle. 
However,  the  technique  discussed  above  whereby  the  surface  fitting  poly¬ 
nomials  are  determined  first  before  any  calculations  are  performed  would 
seem  to  solve  this  problem.  Now  the  interpolations  are  relegated  to  merely 
solving  given  polynomial  expressions  and  the  extra  computational  time  involved 
is  relatively  small.  A  few  less  obvious  benefits  are  gained  by  using  this 
technique.  In  steady  flow,  the  streamline  patterns  are  easily  traced -both 
along  the  body  and  throughout  the  flow  field.  The  output  data  is  also  more 
easily  handled  since  the  flow  properties  are  usually  known  in  planes  orthogonal 
to  the  body  axis,  rather  than  at  random  nondescript  points  throughout  the  flow 
field. 

3.  3  PRACTICAL  NUMERICAL  INTEGRATION  PROCEDURE 

The  preceding  considerations  illustrate  to  some  extent  the 
network  to  be  used  for  calculating  a  point  within  a  flow  field.  However,  in 
order  to  calculate  an  entire  flow  problem  certain  boundary  points  must  be 
included.  In  general,  for  the  class  of  problems  being  discussed  here, 
points  that  lie  on  a  body  surface  and  shock  surface  will  be  required.  In 
the  discussion  that  follows  the  essentials  of  the  integration  process  will  be 
included  for  each  of  these  three  basic  unit  processes.  The  relationship 
which  exists  between  the  three  different  types  of  points  and  the  initial 
surface  is  shown  in  Figure  15. 
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3.3.  1  Starting  Procedure 

The  mesh  points  on  the  initial,  as  well  as  succeeding,  data  planes 
are  systematized  in  terms  of  rays  extending  out  from  the  body.  Each  ray  on 
a  plane  at  X  =  X  is  labeled  O  and  each  point  on  the  ray  l  ,  so  a 

4 

typical  point  is  given  by  .  Points  Pfn  correspond  to  body  points  and 

points  correspond  to  shock  points  where  M  -  max  i  on  ray  o 

The  *.t/  -coordinate  system  will  generally  be  chosen  such  that  the  body 

axi.-  coincides  with  the  X  -axis.  A  new  plane  X-  X  i- A.  is  located  first. 
A  X  is  determined  by  the  C-F-L  stability  condition  which  requires  that 


±1  <1 

A  q  Co 


(72) 


where  Ag  is  the  maximum  distance  allowed  along  a..y  streamline,  Aq 
is  usually  the  minimum  mesh  spacing  and  ^  and  cu  are  the  minimum 
velocity  and  maximum  speed  of  sound  respectively.  It  should  be  noted 
that  a  searching  procedure  is  generally  required  among  the  points  in  the 

data  plane  in  order  that  Equation  (72)  be  satisfied.  However,  this  search 
can  generally  be  confined  to  regions  of  the  flow  field  which  are  a  priori 
known  to  contain  relative  maxima  and  minima  of  the  desired  flow  variables 
(planes  of  symmetry,  etc.).  Given  a  point  in  the  data  plane,  Equation  (72) 
insures  that  the  Mach  forecone  from  a  point  P^  on  the  streamline  through 
P^  is  contained  within  the  eight  initial  data  points  immediately  surrounding 
Pj.  (See  Figure  14)  The  nine  points  shown  determine  the  region  in  the 
initial  plane  over  which  interpolation  surfaces  are  fitted.  If  a  larger  step 
size  A*  is  allowed,  a  correspondingly  larger  region  in  the  initial  data 
surface  must  be  fitted  (with  the  resulting  larger  truncation  errors). 


3.3.2  Field  Point  Routine 

The  solution  at  a  field  point  involves  essentially  six  separate 
processes:  (1)  the  surface  fits  to  the  flow  variables  in  a  region  surrounding 
the  initial  point  are  determined,  (2)  a  new  field  point  is  located, 

(3)  bicharacteristic  lines  leading  from  the  new  point  back  to  the  initial  plane 
are  located,  (4)  the  flow  quantities  at  these  base  points  are  obtained  using 
the  fits  developed  in  step  1,  (5)  the  flow  properties  at  the  new  point  are 

determined  to  a  first  approximation,  and  (6)  the  flow  properties  are  iterated 
to  convergence  to  second-order  in  the  step  size. 
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Step  1 

Given  an  initial  Field  point  (Figure  14)  nine  points  including 
Pj  are  fit  with  an  exact  interpolating  polynomial  for  each  function 
required  in  the  compatibility  equation  written  along  the  bicharacteristics. 
For  the  function  G t-  it  will  probably  be  most  convenient  to  use  the  pressure 
pt,  flow  angles  &  and  ^  ,  Mach  angle  P  ,  the  function  'f  and  the  product 

.  In  this  way  computation  time  can  be  reduced  since  the  compat¬ 
ibility  equation  involves  only  these  functions.  Compared  to  the  concentration 
/Yl  ,  each  function  &■  varies  rather  slowly  throughout  the  shock  layer. 

As  mentioned  above,  this  fact  represents  a  distinct  advantage  of  the 
proposed  scheme. 

The  coefficients  A1 A2t  ,  A3 ■  ,  etc.  to  be  used  in  the  inter¬ 
polation  for  the  function  G,  are  determined  by  solving  the  set  of  equations 

(Gt)j  =  a1l  +  A2c  yj  *  A31?J  +  56* +  A\%  h  +  A6,ch 

(73) 

+  A  7t  ft  y*  +  A  S  jfy  +  A  9-  1  j  *  1, 2,  •  •  •,  9 

written  for  the  nine  initial  points  Pj(z,  yj ,  fj)  . 

Step  2 

The  new  field  point  P^  is  located  as  the  intersection  point  of 
the  streamline  from  Pj  with  the  given  plane  £  =  It  +  A  t>  as 


ft  =  & 


+  Jen  QcLS 


h  *  h  -  c°3  &, <*s 
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where 


cCS 


A* 

cos  B1  cos  <p1 


(77) 


is  the  distance  along  the  streamline. 

If  the  new  plane  is  not  perpendicular  to  the  P  -axis,  Equations 
(74)-(77)  along  with  the  equation  of  the  plane  locate  the  coordinates  of  the 
new  field  point. 

Step  3 

With  the  new  field  point  located,  the  first  approximation  for 
the  flow  properties  there  are  required.  Since  no  calculated  properties  are 

I 

available  at  P  ,  they  must  be  estimated  initially.  A  pressure  gradient 
K  -  _JL — L  is  assumed  at  while  the  remaining  flow  properties  at  P^ 
are  taken  equal  to  their  corresponding  values  at  the  base  point  P^. 

The  compatibility  equation  (71)  corresponding  to  the  Mach 
forecone  from  P^  involves  only  three  dependent  variables  (  f>  ,  6  ,  ip  ). 

It  would  be  reasonable  to  expect  that  only  three  bicharacteristic  lines 
through  P^,  would  be  necessary  to  solve  for  these  three  unknowns. 

However,  while  at  least  three  bicharacteristics  are  indeed  necessary,  it 
has  been  found  that  more  are  required  in  order  to  give  sufficient  accuracy 
to  the  flow  quantities  at  P^.  The  reason  for  this, as  will  be  discussed  in 
more  detail  below,  is  that  in  writing  Equation  (71)  in  finite  difference 
form  along  the  bicharacteristics, certain  higher  order  terms  are  neglected. 
Since  the  iteration  scheme  to  be  described  in  Step  5  below  yields  values 
accurate  to  second  order  in  the  step  size,  in  general,  these  terms  are 
not  negligible. 

Four  bicharacteristics  through  the  new  point  P^  are  chosen. 

The  base  points  P^,  P^,  P,_,  P^  in  the  initial  plane  are  located  by  solving 

Equation  (54)  for  four  values  of  the  parametric  angle  £  -  Ji  J)-  , 

TT  and  .  Typically  when  the  initial  plane  is  given  by  -  Xk,  point  P  , 

c  o 

(  o  =  3,  4,  5,  6)  is  located  from  the  following  equations, 

X,  ^  Xk  (78) 

ij-  »  tjt-(scr)  62  cos /3Z  +  Sin/Sj  cos  Qt  cos  Sc  )  dLt  (79) 
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(80) 


fa  -  }z  -  (cos/S£cos  &zsin</>2  -  sinjd2  \sinty3sin&2c  ?3<S£  -cos  p2sin  B^dL^ 
where 

3?  ,  -  X  • 

<*/.  a - - - — i  a(81) 

‘  cos  &z  cos  <fl2  -  sin fi2  [stn  &2  cos  <f2  cos  6k  +  sin  ip2  s in  6t  J> 

is  the  distance  along  the  respective  bicharacteristic. 

Step  4 

The  next  step  is  to  determine  the  flow  properties  at  the  base 
points  of  the  bicharacteristics.  The  flow  properties  at  each  base  point  can 
be  determined  by  substituting  its  coordinates  into  the  interpolation  formulae 
derived  in  Step  1.  A  new  estimate  of  the  parametric  angle  dL  is  also 
required  at  the  base  of  each  bicharacteristic.  Equations  (78)  -  (81)  above 
with  coefficients  evaluated  at  each  base  point  respectively,  furnish  this  new 
estimate.  The  actual  St-  used  in  the  computations  along  each  bicharacteristic 
then  represents  an  average  between  that  at  point  and  the  associated 
base  point. 

Step  5 

The  compatibility  equations  (56),  (57),  (58),  (59)  and  (71)  are 
written  in  finite  difference  form,  and  used  to  find  first  estimates  to  the 
flow  quantities  at  Point  P  .  Equation  (71)  is  put  in  finite  difference  form 
by  substituting  the  following  approximate  expressions  for  the  derivatives 
of  a  function  f  in  the  bicharacteristic  direction. 


Hi  ,  fi  -  fi 
3L  AL 


c  -  3,  V,  S,  C  (82) 


where 

At  -  \(*z-xt ■Affc-yJ*  '  (83) 


As  a  result,  the  following  linear  equations  for  the  variables  p  ,  9  and  f 
occur  along  each  bicha  -acte r  istic  P{  P^  i  •  3,  ••• ,  6 
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(84) 


~T7  <PM  -  Pc)  *  COS  S„(ez  ■  ei)  *  COS  ei  at>7  6i 

¥  4 

*  -  si»A (cos6>  coa ei  (HI  - sin 6i (IS),  +  )&Ll 

The  only  unknowns  in  Equation  (84)  are  the  quantities  pt  ,  &2  and  <p2 
at  the  new  point  P^.  One  would  thus  expect  that  Equation  (84)  written 
along  three  bicharacteristics  could  be  solved  simultaneously  for  the 
unknowns  at  P~.  However,  before  the  flow  properties  at  point  P9  can  be 
obtained,  the  derivatives  ——  in  the  non-characteristic  direction  N  must 

*7/V 

be  evaluated  numerically.  The  existence  of  derivatives  in  the  direction 

normal  to  a  bicharacteristic  represents  the  major  difference  between 

two-variable  and  three-variable  method  of  characteristics.  In  two- 

independent  variables,  the  compatibility  equations  written  along  the 

characteristics  are  ordinary  differential  equations  whose  difference 

analogs  can  be  solved  immediately  for  the  flow  quantities  at  the  new  point. 

In  this  sense  the  compatibility  equations  are  weaker  for  three -variable 

problems  than  the  corresponding  characteristic  equations  in  two-variable 

problems.  On  the  other  hand,  as  has  been  pointed  out  previously,  there 

are  an  infinite  set  of  such  relations  available  in  the  three -variable  problem. 
2  1 

Butler  has  utilized  some  of  these  extra  relations  in  order  to  eliminate 
the  derivatives  in  the  non-bicharacteristic  direction  at  the  unknown  point. 
Here,  extra  relations  are  introduced  in  order  to  obtain  more  accurate 
average  flow  properties  at  the  new  point.  The  two  techniques  should  be 
essentially  the  same  inasmuch  as  both  yield  improved  approximations  to 
the  derivatives  in  a  non-bicharacteristic  direction. 

Considering  any  three  base  points,  say  P^>  P^  and  P^,  the 
partial  derivatives  can  be  evaluated  from 

<29  _  SO  Sx  +  SO  dy  +  S9_  c ^ 

S/V  '  Sx  SN  +  dj  ~Sn  +  S$  S/V 

where  the  ,  etc.  are  obtained  from  equations  of  the  form 

Sx 

m  +§7{fzK)  <86) 


and  the - -  etc.  are  obtained  from  the  coordinate  transformation  (68)  as 

<2# 


51 


(87) 


-  -  sin  0  cos  u>  sin  £  -  cos£  sin  d/ 

<?A/  r  T 


JJL 


-sin  6  cos  6 


gjy  *  Sin  9  sm  ip  sm  6  +  COSO  cos  <j) 


(88) 

(89) 


Note  that  in  Equation  (86),  9  is  treated  as  an  unknown  so  that  when  (85)  is 

combined  with  (84)  both  the  right  and  left  hand  side  of  (84)  contain  the 

d  V 

unknowns  ^  and  ^  •  The  procedure  for  obtaining  is  similar. 

With  the  coefficients  thus  evaluated,  the  three  linear  equations  of 

the  form  (84)  written  along  the  bicharacter  istics  P^P^,  P^P^  and  P^P^ 

can  be  solved  simultaneously  to  yield  first  approximations  to  the  flow 

variables  p  ,  9  and  V'  (say  P2<0>  >  pf*  )  at  point  P.,.  Three  more 

bicharacteristics  (P^P^,  P^P^  anc*  ^2^6^  are  t^len  c8osen  and  using  the 

technique  described  above  another  estimate  for  the  flow  quantities  p  ,  9 

and  $  at  P  (  p2U>,  9zl) ,  P2U)  )  is  determined.  Finally,  the  first 
c* 

approximation  to  the  pressure  at  P-,  is  given  by 


(90) 


Similarly,  the  average  values  of  0  and  <pz  are  formed.  It  should  be 
noted  here  that  more  sets  of  bicharacteristics  can  be  located  and  still  more 
first  estimates  to  the  flow  properties  at  P^  determined.  Then  Equation  (90) 
would  represent  an  average  over  all  such  estimates.  In  general,  however, 
if  four  bicharacteristics  are  positioned  evenly  around  the  Mach  forecone 
from  P^.the  addition  of  any  other  estimates  at  P^  would  probably  exceed 
the  overall  accuracy  of  the  remainder  of  the  scheme  and  the  final  results 
would  be  changed  only  slightly. 


The  compatibility  equations  (58)  -  (59)  along  the  streamline  PjP^ 
are  integrated  to  yield  the  remaining  flow  quantities  at  point  P^.  In  general 
for  nonequilibrium  flow,  the  temperature  "  ,  concentrations/^  and 


vibrational  energies 


vary  much  more  rapidly  from  point  Pj 
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to  point  than  do  the  pressure  p,  and  flow  angles  8  and  .  Hence, 

the  equations  along  the  streamline  cannot  usually  be  replaced  by  their 

first  differences  and  integrated  directly  as  was  done  above  for  the 
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bicharacteristic  equations.  Following  Wood  et.  al.  and  Curtis  ,  the 
pressure  is  assumed  to  vary  linearly  along  the  streamline 
and  a  fourth-order  Runge-Kutta  technique  is  used  to  integrate  the  given 
coupled  ordinary  differential  equations  in  several  steps  along  P^P£. 
Although  the  Runge-Kutta  method  is  generally  satisfactory,  for  many  cas  s 
of  high-temperature  air  flows  (e.  g.  conditions  where  both  fast  and  slow 
chemical  rates  are  present  in  the  same  region)  it  fails  to  prevent 
instabilities  from  arising.  Treanor  has  recently  developed  an  algorithm 
designed  especially  for  handling  these  so-called  "stiff"  equations  which 
often  occur  in  reacting  and  relaxing  flows.  The  method  which  is  similar 
to  the  Runge-Kutta  technique  includes  terms  which  are  fifth  order  in  the 
integration  step  size.  The  details  are  given  in  Reference  47  and  will  not  be 
repeated  here.  However,  it  should  be  emphasized  that  the  calculation 
time  per  integration  step  is  generally  quite  long  because  several  rather 
complicated  derivatives  must  be  evaluated. 

Step  6 

Having  determined  a  first  approximation  to  the  flow  properties 

at  a  new  field  point  the  last  step  is  to  utilize  an  iteration  process  which 

determines  the  position  of  the  new  point  and  the  flow  properties  there. 

The  iteration  process  employed  here  is  similar  to  the  "modified  Euler 
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method"  or  "Huen's  first  method"  whereby  steps  2  through  5  are 
repeated  using  average  values  of  the  flow'  quantities.  Along  the  streamline 
the  coefficients  involving  8f  and  are  replaced  by  average  values  of  the 
following  form 


Sen  8f  g  )  -  f. sm9i  +Jtn6^n  (91) 

wherever  they  occur  in  Equations  (74),  (75),  (76)  and  (77).  Here,  the 
superscript  refers  to  the  current  iteration  cycle.  Similarly,  in  Equation 
(78)  through  (81)  the  coefficients  defining  the  bicharacteristic  directions  are 
supplanted  by  appropriate  average  values.  As  an  example,  Equation  (79) 
now  becomes 
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(92) 


(sin  &jn  cos /ijn  **+■  smjQz<n  1^coiQ^n  cos  Sa 


■r  sinQ^n  r)cas  ■*•  sin ftl°  ,)cos6^n  cos  S ^in  r^)d^t<n^ 

U  ~  3 +.$<.) 

The  coefficients  of  {  p2  -  pt  )i  (  -  dL  ),  (  ip2  -  <p^ )  and  ALt  in  the  compatib¬ 

ility  Equation  (84)  are  also  replaced  by  the  average  of  their  values  at 
and  ^‘^respectively. 

The  iteration  cycles  may  be  continued  until  pf'n)  -  p~n  °  . 

9^n)  -  &z(n~f)  ,  <P2<n>  ~  <p2<n~r*  ,  etc.  within  the  prescribed  limits. 

Usually  it  is  only  necessary  to  test  the  convergence  on  the  pressure  p  , 
since  experience  here  and  in  the  two  variable  problems  have  indicated  that 
whenever  p  has  converged  the  remaining  flow  quantities  have  converged 
also. 

Using  the  modified  Euler  integration  scheme  described  above 
along  with  the  averaging  technique  presented  in  Step  4,  the  truncation 
error  is  third-order  in  the  step  size;  that  is,  the  procedure  is  accurate  to 
second  order  in  the  step  size  as  is  the  case  with  two-variable  character¬ 
istics,  When  the  averaging  discussed  in  Step  4  is  not  included,  the  process 

is  probably  limited  to  first  -order  in  the  step  size  through  the  assumption 

,  ,  ,  38  ,  r-  ioc\ 

that  the  derivatives  ,  -3 —  and  — —  as  given  by  Equation  \oo)  are 

ax-  <7 

constant  over  the  network.  Thus,  the  addition  of  an  extra  bicharacteristic 
relation  is  seen  to  be  essential  in  order  to  achieve  the  desired  accuracy . 
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3.  3.  3 


Bodv  Point  Routine 


Calculation  of  a  point  on  the  body  surface  involves  essentially 
the  same  six  steps  which  were  necessary  for  a  successful  field  point 
routine.  The  body  surface  will  be  assumed  to  be  given  by 

6  (x,y,2)  =  0  (93) 

B  represents  either  a  known  function  or  a  surface  fitting  element  in  the 
region  of  the  body  point  being  calculated.  In  this  way,  rather  general  body 
shapes  and  motions  can  be  accounted  for.  The  geometry  of  a  typical  body 
point  calculation  is  shown  in  Figure  1  6, 

Step  1 

The  first  step  of  the  body  point  routine  is  exactly  the  same  as 
that  for  the  field  point  procedure.  Nine  initial  points  including  eight  nearby 
the  initial  body  point  P^,  are  fit  with  interpolating  surface  elements  of  the 
form  of  Equation  (73)  for  each  of  the  variables  p  ,  8  ,  p  , ^3  ,  f 

and  .  Note, from  Figure  16, that  of  the  eight  points  chosen, only  two  are 

body  points,  the  remaining  six  being  field  points. 

Step  2 

The  coordinates  of  the  new  body  point  are  obtained  as  the  inter¬ 
section  of  a  plane  through  P^  formed  by  the  body  unit  normal  (  ^  ^  >  n3  ) 

and  the  unit  tangent  velocity  ( )  vector  at  P^  with  the  body  surface 
(1)  and  the  given  plane  X^  =  ^  *  A  * .  Since  the  normal  plane  is 

,  _  _ .  I  1,  n,  ”3  |  - 

(A/  x  S )  '  d  L  =  \Cos  9  cos  ip  s/n  9  Cas  &  Sm  Ip  \  •  c/L 

l  t  X  !  (94) 

This  leads  to  simultaneous  solution  of  the  equations 


0 


FIELD  POINTS 
•  -  BODY  POINTS 
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J 


Figure  16  BODY  POINT  NETWORK  -  THREE  INDEPENDENT  VARIABLES 
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and 


[C/7, ),  c«  0,  cos  ft  -  (n,)t  cos  & sen  ft] 

+  [?/!,),  .s^0,  -(n2)1  cose,  COS  ft]  (}  •},)  {9b* 

-  +  [(»,  ) -(nt)t  coi&i  Sin  ft  ]  (*g  '  *i) 

The  solution  defines  the  body  point  P0  (  ).  Equations  (95)  and  (96) 

can  be  solved  by  using  the  second  order  Newton-Raphson  algorithm  with 
y<tJ-  yt  ,  and  f, ,  serving  as  initial  guesses  to  the  coordinates  of 
point  P^. 

Step  3 

As  in  the  case  of  the  field  point  network,  most  flow  quantities 
at  the  new  body  point  are  initially  taken  equal  to  their  respective  values 
at  point  Pj.  The  initial  guess  for  pressure  will  again  be  given  by  Pz  =  Pf 
(Kt! )  where  X  is  a  known  gradient. 

Bicharacteristic  lines  leading  from  point  back  to  the  initial 
plane  are  determined  next.  Three  bicharacteristics  will  generally  be 
required  here  in  order  to  obtain  sufficient  accuracy  in  the  calculations  at 
P^.  Two  bicharacteristic  elements  are  located  in  the  tangent  plane  to  the 
body  surface  at  P^  and  the  third  in  the  normal  plane  containing  the  stream¬ 
line  P1P2*  E9uations  (78)  -  (81)  serve  to  define  the  base  points 
P  (  i  =  3,  4,  5)  when  S:  at  point  P,  is  known.  Since  the  streamline  P.P_ 
lies  on  the  body  surface,  the  p?rametric  angle  ,  defining  the  bicharacter¬ 
istic  in  the  normal  direction, is  given  by  the  scalar  product  of  the  unit 
normal  vector  (  ,  n3)  and  the  coordinate  vector  rj  (Equation  90) 

shown  in  Figure  1. 

Sj  -  cos~/(~(nl)1  Stn0,  cosftt  (n2  )  cos  0,  -(na),  sts>e,  sen  ft]  (97) 
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m 


Then, 


s*.s  -  s* f  f  ,98) 

define  the  remaining  two  bicharacteristic  directions  through  the  new  body 
point  P.,. 

Step  4 

The  flow  quantities  at  the  base  points  are  determined  from  the 
interpolation  surface  fits  found  in  Step  1.  This  step  is  exactly  the  same 
as  Step  4  in  the  field  point  procedure. 

Step  5 

The  next  step  is  the  determination  of  the  flow  properties  at 

point  P^.  For  the  body  point,  this  involves  the  solution  of  two  compatibility 

equations  together  with  the  boundary  condition  of  flow  tangency  at  P  . 

Choosing  two  of  the  three  bicharacteristics  located  in  Step  3  (say  P  P  and 

3  j  £ 

P^P^),  the  derivatives  ^7  which  appear  in  the  compatibility  equations 

along  each  bicharacteristic  must  be  computed.  These  derivatives  are 

evaluated  in  much  the  same  manner  as  that  used  for  the  field  point  solution. 

Now,  however,  in  evaluating  etc.  in  Equation  (86)  the  initial  body- 

point  Pj  is  used  as  one  of  the  base  points. 

The  compatibility  equation  (84)  along  anc^  are 

two  equations  in  three  unknowns  pt  ,  &2  and  P2  .  The  pressure  p  is 
eliminated  from  each  to  yield  one  equation  relating  &2  and  P2  .  The 
conditions  of  flow  tangency  provides  the  additional  relation  between  and 

^2  as 

(n,)2  cos  &t  cos  p3  +  sen  &2  *  (,ns)2  cos  sin  p2  -  O  (99) 

Using  a  bivariate  Newton-Raphson  procedure  these  two  equations  can  be 
solved  for  initial  estimates  tc  the  flow  angles  at  P^.  The  pressure 
is  then  found  from  one  of  the  original  compatibility  equations.  Repeating 
the  above  for  bicharacteristics  P^P^  anc*  ^5^2  anot^er  estimate  for 
0  ,  ip  and  p  at  P^  is  obtained.  The  final  first  approximation  to  the 
pressure  and  flow  angles  at  P^  is  represented  by  the  average  of  these  two 
solutions. 
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The  remaining  flow  properties  at  are  obtained  by  integrating 
the  compatibility  equations  (57)  through  (59)  along  the  streamline  in  the 
same  manner  as  was  used  in  the  field  point  routine. 

Step  6 

Using  the  "modified  Euler  process"  described  in  Step  6  of  the 
field  point  procedure,  Steps  2  through  5  are  iterated  to  convergence.  The 
location  and  properties  of  point  are  then  stored  and  the  computation 
continued  to  the  next  point.  This  completes  the  body  point  solution  for 
steady  flow. 
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3.  3.  4 


Shock  Point  Routine 


The  shock  point  process  is  the  most  complex  of  the  three  basic 
unit  processes  required,  although  it  follows  the  same  basic  pattern  as  has 
already  been  established.  Now,  however,  not  only  must  the  new  shock 
point  be  located  and  the  composition  and  flow  properties  there  calculated, 
but  the  orientation  of  the  shock  surface  at  the  new  point  must  also  be 
determined.  The  geometry  involved  in  the  shock  point  solution  is  shown 
in  Figure  1  7 . 


Step  1 

This  step  is  exactly  the  same  as  Step  1  of  the  body  point  proced 
ure  with  the  body  points  being  replaced  by  three  shock  points.  The  form 
of  surface  fit  chosen  and  the  flow  variables  fitted  are  the  same  as  those 
given  in  Step  1  of  the  field  point  procedure. 


Step  2 


Since,  unlike  the  body,  the  location  of  the  shock  is  not  known, 
a  point  on  the  shock  will  be  located  as  the  point  of  intersection  of  a  line 
tangent  to  the  shock  surface  through  the  initial  shock  point  Pj  and  the 
plane  X  -  %t+Ax.  The  tangent  line  lies  in  the  plane  formed  by  the  shock 
normal  vector  A/s  ,  and  the  velocity  vector  ^  .  Its  direction  is  .therefore, 
chosen  as  that  of  the  velocity  component  tangent  to  the  shock  surface  from 
the  point  and  is  given  by 


7»  * 

7'*  \l 


-  r,  i  +  r 


V,  J  *  r>1 10 


(100) 


With  r,  established,  the  position  of  the  new  shock  point  is 


Xz  *  X,  +  A/t 


(101) 


= 


(r*)t 


Ax, 


(102) 


h  =  h 


(rf) 


(rx) 


-  Ax 


(103) 
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Figure  17  SHOCK  POINT  NETWORK  -  THREE  INDEPENDENT  VARIABLES 
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Step  3 

The  flow  conditions  behind  the  shock  at  the  new  shock  point 
are  determined  by  solving  the  shock  wave  equations  (Rankine  - Hugoniot 
equations)  using  given  values  of  the  normal  and  the  "free  stream" 

conditions  at  P,.  The  normal  /V,  is  initially  chosen  as 


(104  ) 


where 

(  Ki- 


K, 

~JTL 


is  a  known  gradient  based  upon  the  two  previous  shock  points. 
The  flow  properties  incident  to  the  shock  o  ,  p  ,  /o  , 

rZ  +  2  +  -2  -t 

( )  are  either  known  or  they  can  be  determined  by  inter¬ 


polation  between  the  closest  known  points  in  the  incident  stream.  {In  this 


way,  both  internal  and  external  shock  wave  surfaces  may  be  handled  by  using 
exactly  the  same  technique.  )  Quantities  ahead  of  and  behind  the  shock  are 


denoted  by  a  plus  or  minus  sign,  respectively.  Letting  the  subscripts  N 
and  T  denote  the  normal  and  tangential  velocity  components,  respectively, 


the  shock  equations  can  be  written  in  the  form 


= 

/>-  f»- 

(105) 

= 

P -  +P- 

(106) 

h+±  9  2 

2 

- 

h- 

(107) 

- 

<2+  si  no" 

(108) 

- 

COS  o~ 

(109) 

fr. 

- 

?r. 

(HO) 

= 

(111) 

- 

(112) 
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where  cr  is  the  shock  wave  angle  . 


In  the  general  case,  these  equations  must  be  solved  numerically 

for  conditions  behind  the  shock.  To  this  end,  given  the  shock  wave  angle 

cr  at  the  point  P-,  the  conditions  across  the  shock  can  be  found  using  the 

L  39 

iterative  method  of  Garr  and  Marrone  .  Here,  the  n*-essure  pt_  is 
initially  assumed  and  equations  (105)  -  (107)  above  are  solved  for  the  remaining 

unknowns  p2_,  )  and  h2_  .  The  process  is  then  iterated  by  varying 

the  pressure  until  the  total  enthalpy  behind  the  shock  matches  that  ahead  of 
the  shock  within  a  specified  relative  error.  The  flow  velocity  components 
are  next  calculated  from  the  known  values  of  Ns ^  and  N_  as 


Then, 


U2- 

--  **+  + 

(113) 

*Z- 

£ 

* 

1 

1 

(114) 

•*2- 

A/Sjr(^/V- 

<j«.) 

(115) 

the  flow  angles  (&2  )_  and  (  are  given  by 

(8*)-  =  =J  (116) 

\  V  U  *  4-  7r  2/ 

(*t) 11171 


Thus,  all  flow  properties  behind  the  shock  wave  are  known  at  the  new  shock 
point  P^. 

Step  4 

Next,  the  locations  of  three  bicharacteristics  leading  from  the 
new  shock  point  P^  back  to  the  initial  plane  X  =Xh  are  desired.  In  a 
similar  manner  to  that  used  in  the  body  point  solution,  two  of  the  bicharacter¬ 
istics  are  selected  to  lie  in  the  shock  surface,  while  the  third  is  chosen  to 
intersect  the  initial  surface  at  a  point  in  the  field.  The  base  points  P^ 
and  P^  on  the  shock  surface  in  the  initial  plane  are  determined  as  the 
points  of  intersection  of  the  Mach  cone  (Equation  (48)  with 
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the  shock  surface  and  the  initial  plane  7C  =  A'A.To  this  end,  the  element 
P^PjPg  (see  Figure  17)  of  the  shock  surface  in  the  initial  plane  is  fitted 
with  a  quadratic  in  y  and  as 


y  ,  Ay*  4-  6<yt  D  (lly) 

Then,  equations  (118)  and  (48)  yield  the  coordinates  of  the  two  base  points 
P^  and  P^.  The  parametric  angles  defining  the  bicharacteristics  P^P^ 
and  P^P^  can  he  obtained  from  Equation  (50)  with  coefficients  evaluated  at 
P_  and 


+  (y, -yi)2+(h -*;)*) Vt  is3^  (n9) 

Then,  a  third  bicharacteristic  is  located  midway  between  P^P^  and  P^P^ 
on  the  Mach  forecone  by  taking 


at  point  P^  and  solving  Equation  (54)  for  the  coordinates  . 

Step  5 

This  step  is  exactly  the  same  as  Step  4  in  the  field  and  body 
point  procedures.  Given  the  coordinates  of  base  points  P7,  P^  and  P,_ 
the  function  surface  fits  evaluated  in  Step  1  can  be  used  tc  determine  the 
base  point  flow  properties. 

Step  6 

In  Step  3,  all  flow  quantities  behind  the  shock  at  point  P^  were 
determined  by  solving  the  Rankine -Hugoniot  equations.  These  flow  proper¬ 
ties  must  also  satisfy  the  compatibility  equations  written  along  the 
bicharacteristics  P^P^i  P^P^  anc^  ^5^2  w^en  the  shock  wave  is  correctly 
located  and  oriented  at  P^.  This  fact  can  be  used  to  develop  an  iteration 
scheme  to  converge  the  position  and  orientation  of  the  shock  wave  element 
through  P^. 
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Since  conditions  behind  the  shock  at  any  point  cannot  be  multi¬ 
valued,  whenever  one  flow  property  is  correct  the  remaining  properties 
must  be  right  also.  Consider  that  the  flow  angles  )_  and  ( liL  )_  as 
determined  in  Step  3  above  ,are  correct.  Then,  each  of  the  three  compat¬ 
ibility  equations  written  along  the  bicharacteristics  P^P^,  P^P£  anc*  ^5^2 
can  be  solved  to  yield  three  values  of  the  pressure  pt<^  at  P^.  The 
partial  derivatives  which  occur  in  the  compatibility  equations  are 

Off 

evaluated  here  in  the  same  manner  as  was  used  in  the  field  and  body  point 
procedures.  (Equations  (85)  -  (89).  ) 

Comparing  the  average  pressure  px  *  as  obtained  from  the 
compatibility  equations  with  the  pressure  p  from  the  shock  relations,  a 
new  estimate  to  the  shock  normal  at  P^  can  be  determined  by  defining 


e'"1  = 


(121) 


Then,  new  values  of  the  direction  cosines  of  the  normal  N.  are 

5i 

selected  by  linear  interpolation  as 


■  « 


(n-t)  tr>) 


(n)  (n-f)\ 


V<' 


'  (n) 

e  -e 


(122) 


’> )  ,12  3) 

while, 

",r* ■  •  V  -  ■ n\  Y/a  1 1 24) 

<fz  y 

Again,  thf*  superscript  ( n )  denotes  the  iteration  cycle.  This  procedure  is 
repeated  i.  an  iteration  process  using  average  flow  properties  and  coeffic¬ 
ients  to  locate  the  new  shock  point  and  the  bicharacteristic  base  points  as 
in  the  field  and  body  networks.  Average  values  are  also  used  for  the 
coefficients  of  the  respective  compatibility  equations  associated  with  the 
bicharacteristics  ^4^2  anc*  ^5^2"  Iterat'on  continues  until 
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*  Test  Constant 


(125) 


ft-  '  ?t- 
Pt- 

at  which  time  the  position  and  orientation  of  the  shock  wave  are  found  that 
yield  flow  properties  satisfying  the  compatibility  equations  to  the  desired 
degree  of  accuracy.  This  completes  the  shock  point  solution. 

3.4  SIMPLIFIED  FLOW  MODELS 

The  methods  derived  above  apply  to  three  dimensional  steady 

gas  flows  in  which  chemical  reactions  and  vibrational  relaxations  occur  at 

35 

finite  rates.  Experience  with  finite  rate  programs  indicates  that  even 

for  calculations  involving  two  independent  variables  the  full  equations  are 

apt  to  require  long  machine  times  -  especially  in  near  equilibrium  flow 

situations.  Therefore,  it  may  not  always  be  desirable  or  even  necessary 

to  consider  the  full  nonequilibrium  flow  equations.  Instead,  following 
49 

Curtis  ,  appropiiate  simplifications  of  the  chemistry  may  be  introduced, 
Specifically,  three  sets  of  simplifying  assumptions  relative  to  the  chemistry 
and  the  thermodynamics  of  the  gas  flow  are  presented  here.  Each  one  is 
listed  below,  along  with  the  procedures  required  for  carrying  it  out. 


3.  4.  1 


Ideal  Gas  Model 


In  this  case,  all  chemical  production  terms  are  equated  to  zero 
and  the  specific  heat  is  constant  corresponding  to  no  excitation  of  the 
vibrational  degrees  of  freedom.  Therefore,  the  variables  and  terms 
involving  them  may  be  eliminated  from  all  preceding  equations:  RS ,  N  V  t 

r— — * 

P' -  ,  £■  ,  /f^r,  ,  Q;  and  A"  .  The  speed  of  sound  is  computed  from 

*/  w  J 

(126) 

and  for  steady  flow  the  velocity  is  given  by 

jpii  -  (H^-CpT)  ( 1 27) 

3.4.2  Frozen  -  Inhomogeneous  Model 

The  frozen-inhomogeneous  model  corresponds  to  the  case  where 
the  gas  composition  along  streamlines  (or  particle  paths)  has  frozen  in  a 
partially  dissociated  state.  The  composition  does  vary  across  streamlines 
however,  and  thus  throughout  the  shock  layer  the  gaseous  mixture  is 
inhomogeneous.  Assuming  vibrational  equilibrium,  frozen  flows  can  be 
calculated  by  setting  the  following  terms  equal  to  zero  in  the  compatibility 
equations:  AJV  ,  R^  ,  and  t  .  The  frozen  specific  heat 

is  obtained  from  Equation  (28)  with  A/f  -  and  the  frozen  speed  of 
sound  is  again  given  by  (34).  The  netword  of  intersections  of  streamlines 
and  reference  planes  is  seen  to  be  suitable  for  calculating  frozen  flows 
because  the  concentrations  of  the  various  species  are  forced  to  remain 
constant  along  each  streamline. 

3.4.3  Equilibrium  Model 

The  equilibrium  model  assumes  that  local  thermochemical 
equilibrium  will  be  established  at  every  point  in  the  flow  field.  Again, 

Rv-j  >  6j  ,  A/V  ,  Q t-  ,  and  t  are  set  equal  to  zero  in  the 

compatibility  equations.  Two  distinct  approaches  are  possible  here. 

First,  the  temperature,  density  and  speed  of  sound  could  oe  determined 
from  polynomial  fits^  or  special  curve  fits^^  of  tabulated  air  data. 
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In  this  case,  the  equilibrium  rather  than  the  froze"  speed  of  sound  is  used. 
Knowing  the  pressure  and  temperature  at  each  point  in  the  field  the 
equilibrium  concentrations  can  be  determined  uniquely.  A  second  technique 
would  be  to  use  an  iterative  method  to  solve  the  full  coupled  equilibrium 
equations  at  every  point  in  a  manner  similar  to  that  used  by  Boyer 
While  the  latter  scheme  is  more  accurate,  the  first  method  would  involve 
far  less  computational  time  and  should  probably  be  used  in  any  multi¬ 
dimensional  characteristics  program.  The  accuracy  of  the  approximate 
method  can  generally  be  held  within  1%  which  is  sufficient  for  almost  all 
practical  purposes. 

3.  5  INITIAL  VALUE  SURFACE 

The  body  point,  field  point  and  shock  point  solutions  described 
above  can  be  used  as  building  blocks  from  which  the  complete  flow  field  can 
be  constructed.  However,  before  the  method  of  characteristics  can  be 
utilized  to  solve  such  flows,  an  initial  value  surface  on  which  the  flow 
properties  are  completely  specified  must  be  provided.  In  many  cases, 
obtaining  these  starting  values  is  not  a  simple  problem.  In  the  discussion 
that  follows  several  possible  procedures  which  provide  the  required 
starting  regions  for  rather  specific  problems  are  discussed. 

In  the  case  of  a  completely  general  three-dimensional  body,  a 

plane  of  inputs  along  which  the  flow  is  initially  supersonic  is  required. 

Unfortunately,  since  the  necessary  data  will  in  general  be  asymmetric  ,  a 

method  for  obtaining  inputs  in  this  general  case  does  not  currently  exist. 

Consequently,  more  specialized  solutions  must  be  attempted.  One  such 

specialization  consists  in  solving  the  flow  field  in  the  vicinity  of  a  spherically- 

capped  body  at  angle -of-attack.  In  this  case,  initial  data  in  a  plane  normal  to 

the  flow  direction  is  axisymmetric  and  may  be  provided  /  using  any  of  the 

well-known  axisymmetric -blunt  body  solutions.  Among  these,  either  the 
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direct  method  using  integral  relations  or  an  iterated  inverse  method 
would  yield  initial  data  which  is  sufficiently  accurate  for  continuing  the 
solution  using  the  method  of  characteristics.  The  only  restriction  to  such  an 
approach  is  that  for  the  angles  of  attack  considered,  the  body  sonic  point  on 
the  compression  side  of  the  body  must  be  located  on  the  spherical  nose. 
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Certain  three-dimensional,  steady  blunt  body  flows  may  also  be 
solved  by  treating  them  as  limiting  solutions  in  time  of  a  properly  posed 
unsteady  flow  problem.  These  methods, which  have  recently  been  advanced 
in  the  literature  ’  '  .consist  first  of  moving  from  a  known  or  attainable 
steady  flow  solution  to  another  desired  steady  solution  using  a  prescribed 
body  motion.  The  body  motion  may  consist  of  either  a  systematic  warping 
from  an  initial  shape  (e.  g,  ,  a  spherical  cap)  to  a  final  shape  (e.  g.  ,  an 
ellipsoidal  nose)  or  of  moving  a  given  body  from  one  configuration  (e.g.  , 
flow  at  zero  angle  of  attack)  to  another  (e.  g.  ,  flow  at  high  angle  of  attack) 
or  any  combinations  of  these  two.  The  final  steady  flow  is  obtained  by 
continuing  the  calculations  in  time  with  the  desired  body  surface  unchanged 
until  the  flow  transients  due  to  the  prescribed  motions  have  disappeared. 

Since  as  has  been  pointed  out  before,  the  governing  equations  for  unsteady 
flow  are  hyperbolic,  the  method  of  characteristics  for  unsteady  flow 
may  be  conveniently  used  as  a  basis  for  solving  the  steady  flow  problem 
using  this  approach.  Here,  starting  with  a  known  flow  field  in  the  initial 
hyperplane  t  -  ta  ,  (e.g.,  the  steady  flow  solution  about  a  three-dimensional 
nose)  a  series  of  unsteady  three-dimensional  flow  problems  are  solved  step- 
by-step  in  time  until  the  flow  field  no  longer  changes  with  time.  Although 
these  computations  may  indeed  be  very  tedious  and  involve  long  computer 
run  times,  the  flow  fields  about  some  rather  arbitrary  bodies  can  be 
determined. 
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Section  4 

RESULTS  AND  CONCLUSIONS 

4.  1  PROBLEM  AREAS  AND  THEIR  SOLUTIONS 

Before  presenting  the  results  obtained  from  the  three-dimensional 
steady  flow  method  of  characteristics  program,  it  would  seem  to  be  desirable 
to  first  discuss  some  of  the  difficulties  which  arose  during  the  check  out  and 
debugging  phases  of  the  program  development.  This  discussion  should  prove 
helpful  to  those  who  attempt  to  utilize  and/or  modify  the  basic  method  as  it 
was  presented  in  Section  3.  The  overall  problems  were  stability,  accuracy, 
and  computer  running  time,  in  that  order  of  importance.  These  problems 
did,  however,  prove  to  be  interrelated,  since  once  greater  accuracy  was 
attained  the  computer  run  time  improved  considerably. 

4.  1 .  1  Interpolation 

As  discussed  previously,  when  the  integrations  are  performed  by 
tiacing  streamlines  step-by-step  through  the  flow  field  the  points  in  each 
plane  normal  to  the  X  -axis  do  not  maintain  a  uniform  mesh  spacing.  Con¬ 
sequently,  what  was  required  was  a  bivariate  interpolation  procedure  which 
was  applicable  to  a  net  of  irregularly  distributed  points.  Unfortunately, 

because  of  the  generality  of  the  problem,  very  little  attention  beyond  that  of 
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jus t  stating  the  difficulties  involved,  ^  '  has  been  given  to  it  in  past  work. 
Generally,  using  polynomials,  only  two  basic  approaches  to  this  problem  are 
deemed  useable.  The  first  possibility  would  be  to  use  a  method  of  least 
squares  approximation  whereby  the  number  of  base  points  through  which  the 
fit  is  to  be  passed  exceeds  the  number  of  coefficients  which  are  required  in 
the  fitting  polynomial.  Using  this  procedure  the  sum  of  the  squares  of  the 
errors  represented  by  the  differences  between  the  value  of  a  given  function  at 
a  point  and  the  value  obtained  for  that  function  from  the  fit  should  be  a  min¬ 
imum.  The  other  possibility  is  to  use  an  exact  polynomial  fit  through  the 
given  base  points.  By  "exact"  here  is  meant  that  the  error  between  the  actual 
value  of  the  function  and  the  values  obtained  from  the  fit  is  zero  at  least  to 
within  the  accuracy  of  the  machine  computations.  The  form  of  each  of  these 
fits  depends  in  turn  upon  the  relative  orientations  of  the  base  points,  the 
number  of  base  points  chosen  and  the  form  of  the  basic  independent  set  of 
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basis  functions  which  comprise  the  fit.  For  example,  for  linear  fit  in  one 
variable,  X  ,  the  basis  functions  would  be  f  and*). 

Both  of  the  above  methods  were  investigated  for  use  with 
the  three-dimensional  characteristics  program.  While  these  studies  were 
were  not  exhaustive,  a  few  definite  conclusions  seemed  to  be  justified. 
First,  for  bodies  with  "rounded"  cross-sections  (i.  e.  ,  axisymmetric  or 
elliptical  cross-sections  which  are  not  too  flat),  as  would  be  expected, 
better  fits  could  be  obtained  by  using  combinations  of  the  polar  coordi¬ 
nates  (rfbJ)  as  the  basis  functions  for  the  fits  rather  than  the  cartesian 
coordinates  (u  ,  y).  This  was  especially  true  in  the  case  of  the  flow 
angles  where  relatively  large  changes  occur  around  the  body.  For  bodies 
which  have  locally  two-dimensional  regions,  it  would  appear  that  probably 
a  better  set  of  basis  functions  would  be  combinations  of  the  cartesian 
coordinates  {X  ,  (j).  This  point  has  not  yet  been  thoroughly  checked  out 
however,  and  probably  even  for  these  areas  the  polar  fits  will  be  satis¬ 
factory.  If  not,  the  program  can  readily  be  changed  to  allow  for  using 
different  basis  functions  over  different  regions  of  the  flow  field.  In  the 
case  of  an  elliptical  cross-section  this  could  best  be  accomplished  by 
translating  the  center  of  the  (  ^  ^  )  coordinate  system  from  the  center  of 
the  ellipse  to  the  focus,  and  then  using  the  (  Y~ ,  CO  )  -  fit  over  the  region  of 

/  )  -  fit 

over  the  remaining  part. 

The  accuracy  of  each  surface  fit  was  improved  considerably  by 
scaling  the  independent  variables  with  respect  to  the  coordinates  of  the 
central  point  used  in  the  fit.  In  many  cases,  improvements  in  accuracy 
amounting  to  two  significant  figures  were  obtained  by  utilizing  this 
technique.  Finally,  lor  the  ranges  of  points  tested  (e.  g.  ,  nets  of  9  and  25 
points  centered  around  a  base  point  in  the  initial  plane)  with  the  polar  fits, 
the  method  using  least  squares  and  that  using  the  exact  fit  yielded  approxi¬ 
mately  the  same  accuracies  for  all  functions  fit.  However,  since  the  exact 
fit  involves  fewer  operations  and  does  in  fact  fit  the  data  at  the  given  points 
exactly,  this  method  was  finally  chosen  for  use  with  the  program.  The 
justification  for  using  an  exact  fit  lies  also  in  the  fact  that  the  data  at  each 
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the  field  covering  the  "rounded"  portion  of  the  ellipse  and  the  , 
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known  point  is  itself  in  a  sense  "exact"  for  a  given  problem  and  hence, 
smoothing  procedures  using  least  squares  or  other  methods  should  not  be 
required.  However,  as  was  indicated  above,  these  results  are  somewhat 
preliminary  and  now  that  the  characteristics  program  is  working,  a  more 
thorough  investigation  of  the  surface  fitting  procedures,  as  applied  to 
several  different  flow  fields,  should  be  undertaken. 

4.1.2  Flow  Angle  Derivatives 

A  more  fundamental  difficulty  occurred  while  attempting  to  apply 
the  characteristics  procedure  to  compute  some  simple  two-dimensional  and 
conical  flow  fields.  Two  flow  problems  were  studied  initially,  (1)  a  simple 
wedge  and  (2)  a  sharp  cone  -  both  at  zero  angle  of  attack.  In  each  case,  the 
minimum  number  of  bicharacteristics  (three)  necessary  to  determine  the 
flow  properties  at  a  new  field  point  was  utilized.  The  wedge  flow  was  success¬ 
fully  reproduced  to  eight  places  in  all  flow  variables  during  each  axial  step. 
This,  however,  was  not  a  very  demanding  test  because  the  flow  field  is 
initially  uniform  and  no  iteration  was  required  in  order  to  determine  the  flow 
properties  at  a  new  point.  Also,  the  normal  derivatives  in  the  compatibility 
equations  were  in  this  case  identically  zero.  The  flow  cvei  a  10°  right 
circular  cone  offered  a  more  stringent  test  on  the  overall  calculation  proce¬ 
dure.  Although  the  initial  data  was  axisymmetric  to  eight  places,  the 
results  in  this  case  were  non-axisymmetric.  The  degree  of  non-axial 
symmetry  varied  according  to  the  gradients  existing  in  the  flow  and  the 
step  size.  However,  the  pressure  and  flow  angles  were  generally  constant 
to  only  three  significant  figures  (i  to  2%)  around  the  cone  after  just  one  or 
two  axial  steps.  These  azimuthal  variations  in  the  flow  properties  around 
an  axi- symmetric  body  could  be  attributed  to  the  fact  that  as  the  calculations 
progressed  around  a  given  ring  o.'  data  points  (  r  =  constant),  the  relative 
orientations  of  the  three  bicvaracteristic s  which  were  used  to  determine  the 
flow  properties  with  respect  tu  the  azimuthal  plane  through  the  initial  data 
point  were  different.  Thus,  the  conditions  at  each  new  point  would  be  expected 
to  be  slightly  different,  in  that  conditions  at  the  respective  base  points,  as 
determined  from  the  interpolations,  were  not  the  same.  The  variations  which 
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did  in  fact  occur  were  surprisingly  large.  These  large  variations  were  attri¬ 
buted  to  the  fact  that  some  first  order  effects  were  included  in  the  compatibility 
equations  through  the  assumption  that  the  derivatives  of  the  flow  angles  with 
respect  to  the  coordinate  directions  were  constant.  Hence,  in  order  ot  obtain 
meaningful  results,  better  approximations  to  the  derivatives  were  required. 
These  were  obtained  by  introducing  extra  bicharacteristics  at  each  point  and 
averaging  the  flow  properties  at  the  new  point  in  the  manner  set  forth  in 
Section  3.  Actually,  in  the  program  three  sets  of  bicharacteristics  (initial 
set  is  rotated  40  s)  are  utilized  to  obtain  this  average.  The  averaging  of  the 
flow  properties  at  the  new  point  improved  the  accuracy  of  the  calculation  by 
as  many  as  two  signifigant  figures.  In  addition,  the  more  accurate  proceedure 
converged  in  three  cycles  rather  than  the  seven  required  without  averaging; 
thus,  although  more  calculations  were  required  within  each  cycle,  the  net  run 
time  was  not  affected  by  the  averaging  procedure. 

4.1.3  Mach  Cone  Intersects  the  Shock 

Another  difficulty  was  encountered  when  a  new  ring  of  field  points 
was  added  near  the  shock  wave  as  a  result  of  introducing  new  streamlines 
from  the  previous  shock  point.  In  this  case,  a  portion  of  the  Mach  forecone 
from  a  new  field  point  intersects  the  shock  wave  surface.  Thus,  the  base 
points  of  some  of  the  bicharacteristics,  which  are  used  to  calculate  the  flow 
properties  at  the  new  field  point,  intersect  the  initial  data  plane  on  the  free- 
stream  side  of  the  shock  surface.  Two  different  methods  were  utilized  in 
attempting  to  solve  this  problem.  First,  the  relevant  portion  of  the  shock 
surface  and  the  flow  quantities  behind  the  shock  were  fitted  by  surface  fits 
similar  to  those  used  in  interpolating  for  the  base  points .  These  fits  were  used 
to  determine  the  coordinates  of  and  flow  properties  at  the  base  points  of  the 
bicharacteristics  which  intersected  the  shock  surface,  while  those  bicharacteristics 
on  the  portion  of  the  Mach  cone  whichdid  not  intersect  the  shock  were  handled  in 
the  usual  manner.  Even  though  this  procedure  seemed  to  work  proDerly,  the 
results  obtained  were  not  compatible  with  the  accuracy  of  the  solution  in  the 
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remainder  of  the  flow  field.  For  an  axisymmetric  flow  field,  the  results 
along  the  new  ring  were  axisymmetric  to  only  three  places,  while  the  rest 
of  the  flow  field  remained  so  to  six  significant  figures.  The  reason  for  the 
discrepancy  is  thought  to  be  due  to  the  fact  that  the  base  points  associated 
with  each  of  the  two  separate  fitting  schemes  are  treated  with  unequal 
weights  since  not  the  same  number  of  bicharacteristics  intersect  the  shock 
as  intersect  the  field.  Time  limitations  prevented  further  development  of 
this  method.  Instead,  another  much  simpler  technique  was  employed.  In 
this  scheme,  the  C-F-L  stability  criterion  was  in  principle  violated  by 
allowing  some  of  the  bicharacteristics  to  pass  through  the  shock  surface 
and  intersect  the  initial  plane  outside  of  the  shock  layer.  The  flow 
properties  at  these  base  points  are  determined  by  extrapolation,  utilizing 
the  interpolating  surface  fits  which  had  already  been  set  up  surrounding  the 
initial  shock  point.  This  procedure  yielded  results  which  were  compatible 
with  the  remainder  of  the  flow  field  solution,  and  is  the  scheme  that  is 
presently  being  used  with  the  program.  Note,  that  violating  the  stability 
criterion  in  a  small  region  does  not  seem  to  cause  course  instabilities  to 
arise.  This  is  because  basically,  instability  is  the  result  of  a  cumulative 
growth  of  errors  throughout  a  calculation  and  does  not  usually  result  from 
a  few  isolated  errors,  as  is  the  use  here.  In  order  to  prevent  the  extra¬ 
polations  from  being  extended  too  far  outside  the  shock  wave,  the  step  size 
is  halved  before  any  points  in  the  new  plane  are  determined,  whenever  a 
new  ring  of  field  points  is  to  be  added.  Note  that  since  the  streamline  is 
inclined  toward  the  shock  wave,  only  a  small  portion  of  the  Mach  forecone 
from  a  new  point  on  the  streamline  will  intersect  the  shock  when  the  step 
size  ,  is  small.  Thus  the  extrapolation  desc ribed  should  be  reasonably 
accurate  as  the  results  obtained  thus  far  using  this  technique  seem  to 
indicate. 

4.1.4  Problems  Associated  with  Large  Angles  of  Attack 

In  general,  as  discussed  above,  orienting  the  X  -axis  (the  axis 
along  which  the  calculations  proceed)  along  the  free  stream  wind  direction 
has  the  advantage  that  the  initial  input  data  is  axisymmetric  and  can  be 
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obtained  from  two-independent  variable  programs.  However,  at  large  angles 
of  attack  (say  30°)  the  flow  in  the  shock  layer  on  the  windward  side  is  at  a 
relatively  low  Mach  number.  Consequently,  the  Mach  lines  are  steeply 
inclined  and  the  step  size  Ax  between  successive  data  planes  is  small.  In 
fact,  whenever  the  sum  of  the  local  two-dimensional  flow  angle  ^  and  the 
Mach  angle  AX  is  greater  than  90°  (see  Figure  18)  the  step  size  AX  becomes 
negative  and  the  calculations  would  proceed  upstream.  The  program  is  pre¬ 
sently  incapable  of  handling  this  situation.  It  should  be  emphasized  however, 
that  this  limitation  does  not  represent  a  deficiency  in  the  overall  method,  but 
is  rather  a  problem  due  to  the  necessity  of  aligning  the  x  -direction  with  the 
free  stream-velocity  vector  in  order  to  utilize  axisymmetric  inputs.  When¬ 
ever  a  method  of  obtaining  non-axially  symmetric  inputs  becomes  available 
the  orientation  of  the  X  -axis  can  be  changed  so  that  the  condition  <3  +  X  =  90° 
is  less  likely  to  occur. 

When  the  axial  step  is  positive  but  small  due  to  the  low  Mach 
number  of  the  flow  in  the  shock  layer,  calculations  of  complete  flow  fields 
would  be  speeded  up  considerably  by  reorienting  the  x  -axis  parallel  to  either 
the  body  axis  or  to  the  windward  side  of  the  body.  This  idea  was  investigated 
by  rotating  successive  data  planes  by  a  small  angular  increment  A/a  about  a 
point  on  the  windward  shock  (Figure  18),  followed  by  a  coordinate  transforma¬ 
tion  to  a  body-axis  system.  Unfortunately,  programming  difficulties  and 
lack  of  time  forestalled  the  realization  of  this  scheme, 

4.  2  TYPICAL  OPERATIONAL  EXPERIENCE  WITH  THE  PROGRAM 

The  results  of  typical  field,  body  and  shock  point  calculations 
are  given  in  Appendix  B,  Tables  B-I,  B-II,  and  B-III.  The  number  of 
iterations  shown  are  typical  for  calculations  to  date  with  the  field  and  body 
generally  requiring  three  to  five  cycles  depending  upon  the  gradients  and  the 
axial  step  size,  while  the  shock  requires  at  least  five.  Convergence  at  a 
shock  point  is  slower,  because  the  first  two  cycles  involve  first,  an  estimate 
of  the  direction  of  the  shock  normal,  and  then  a  given  small  change  in  this 
direction  in  order  to  numerically  determine  the  derivatives  in  the  Newton-Raphson 
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Figure  19  ILLUSTRATION  OF  DIFFICULTIES  ENCOUNTERED  IF  FLOW  MACH  NUMBER  IS  TOO  LOW 


procedure  used  to  obtain  the  next  values  of  shock  parameters.  Thus,  the 
actual  internal  iteration  process  does  not  begin  at  the  shock  point  until  the 
third  cycle  of  the  overall  iteration  has  been  reached.  It  then  usually  takes 
three  more  cycles  to  complete  the  shock  point  computation.  The  results 
indicate  that  the  coordinate*  and  flow  angles  generally  converge  faster  than 
do  the  pressure  and  velocity,  which  usually  converge  at  the  same  rate. 

Actually,  just  the  pressure  has  been  used  for  the  convergence  tests  thus  far 
in  the  program  and  from  the  results  obtained  this  seems  to  work  quite  well; 
however,  in  some  cases,  it  might  be  better  to  test  on  both  the  velocity 
magnitude  and  the  pressure  -  especially  if  a  smaller  convergence  criterion 
than  10”3  is  used.  It  should  be  noted  that  although  the  results  in  Tables  B-I 
through  B-III  are  given  in  terms  of  eight  digit  numbers  corresponding  to  the 
number  of  digits  the  machine  carries,  only  the  first  four  or  five  digits  are 
significant,  due  to  the  various  truncations  which  occur  during  the  computation. 

The  total  computation  time  required  to  completely  solve  for  the 
coordinates  and  flow  properties  at  a  particular  field,  body  or  shock  point  is 
difficult  to  determine  accurately.  It  depends  upon  the  gradients  which  exist 
in  the  flow,  the  spacing  between  data  points,  the  axial  step  size  and  the 
convergence  test  criterion  utilized.  An  estimation  of  the  time  required  can 
be  obtained  by  dividing  the  total  time  taken  to  calculate  several  planes  by  the 
total  number  of  points  calculated.  Although  this  time  includes  the  time  required 
to  read  input  data  from  the  tape,  set  up  the  initial  data  surface,  and  write 
data  back  on  output  tape  as  desired,  the  significant  portion  consists  of  the 
actual  computation  time  -  especially  if  several  planes  (1000  to  2000  points) 
are  computed.  Using  this  technique,  the  time  ranged  from  0.  30  to  0.  50  seconds 
per  point  on  the  IBM  7094  computer.  The  comparable  time  for  the  IBM  7044 
would  be  about  twice  this  value.  Thus,  it  is  absolutely  necessary  to  decrease 
this  time  significantly,  if  many  problems  are  going  to  be  run  using  the  present 
method.  This  will  be  especially  true  when  the  method  is  extended  to  include 
nonequilibrium  thermochemical  effects. 
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RESULTS  FOR  COMPLETE  FLOW  FIELDS 


4.  3 

Three  different  afterbody  flow  cases  were  computed  using  the 
three-dimensional  method  of  characteristics  program.  These  consisted  of 
the  flow  fields  associated  with  a  spherically  blunted  15°  cone  at  both  0°  and 
10°  angles -of-attack,  and  a  spherically  capped  elliptical  afterbody  whose 
cross-sections  were  ellipses  with  the  eccentricity  of  each  varying  as  a 
function  of  the  axial  distance  X  at  0°  angle -of-attack.  Since  the  flows  around 
the  cone  at  0°  angle-of-attack  and  the  elliptical  body  are  symmetric  with 
respect  to  the  y  and  ^  axes,  it  is  sufficient  in  these  cases  to  calculate  only 
one  quadrant  of  the  entire  field.  In  the  case  of  the  blunted  cone  at  angle-of- 
attack  however,  there  exists  only  one  plane  of  symmetry  (the  plane  containing 
the  free  stream  velocity  vector  and  the  y -axis)  and  a  full  180°  of  the  field 
must  be  computed. 

All  computations  were  performed  using  the  same  input  conditions 
obtained  in  the  manner  discussed  in  Paragraph  4.  2.  The  radius  r0  ,  of  the 
nose  sphere  in  each  case  was  0.28  cm  with  the  center  of  the  reference  (xty,^) 
coordinate  system  located  at  the  center  of  the  sphere.  An  initial  line  of 
nine  data  points,  equally  spaced  between  the  body  and  the  shock,  was  provided 
by  the  two  independent  variable  characteristics  program  at  a  non-dimension- 
alized  (with  respect  to  )  axial  coordinate  x  equal  to  -0.4534.  This  initial 
data  line  was  parallel  to  the  y-axis  in  the  -plane.  To  provide  an  initial 
plane  of  data  points  the  initial  line  was  rotated  about  the  T-axis  in  given 
increments  J  uJ  of  the  pohr  angle /i/.  A  different  number  of  meridional  rays 

were  used  for  each  of  the  problems  studied.  Ihe  free-stream  conditions 

.  7 

were  in  all  cases:  -  8.  33,  p#  =3.153x10  ,  -  250 °K  corresponding 

to  an  altitude  of  200,000  feet.  Results  for  each  case  studied  will  be  presented 
separately. 

4.  3.  1  Blunted  Cone  At  Zero  Angle -of- Attack 

Calculation  of  the  flow  field  about  a  blunted  cone  at  zero  angle-of- 
attack  was  carried  out  using  five  initial  meridional  rays  located  22.  5°  apart. 

It  took  approximately  205  steps  to  reach  a  distance  of  8  nose  radii  down  the 
body.  The  computer  run  time  was  nearly  90  minutes  on  the  IBM  7094  computer. 
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While  this  run  time  is  long,  it  is  far  from  the  optimum  that  could  be  obtained 
using  the  program  even  in  its  present  state  of  development.  Since  limited 
computer  time  was  available  for  each  of  the  flow  calculations  described  here, 
it  was  decided  to  allow  the  machine  to  choose  the  axial  step  size  based  upon 
the  minimum  mesh  spacing  in  the  field  so  that  stability  would  be  assured, 
and  the  program  would  proceed  as  far  as  possible  without  the  necessity  of 
re-running  a  problem.  While  this  objective  was  accomplished,  the  step 
sizes  the  program  chose  were  generally  much  smaller  than  would  be  necessary 
for  stability.  This  was  because  the  minimum  mesh  spacing  generally  occurs 
between  the  last  point  on  the  ray  and  the  point  on  the  shock  surface.  When 
new  streamlines  are  added  at  the  shock,  this  mesh  size  is  generally  an  order 
of  magnitude  smaller  than  the  spacings  between  the  rest  of  the  points  in  the 
field.  Consequently,  the  step  size  used  after  adding  points  near  the  shock  is 
an  order  of  magnitude  smaller  than  the  maximum  that  would  probably  be 
used  in  such  computation.  It  usually  takes  five  or  six  axial  steps  for  the  step 
size  to  again  approach  a  reasonable  value.  A  possible  solution  to  this  would 
be  to  set  the  axial  step  size  at  a  given  value  based  upon  the  C-F-L  condition, 
and  then  calculate  several  steps  using  this  constant  value.  Then  the  step 
size  could  be  changed  again  and  another  few  steps  could  be  computed,  and 
so  on.  If  the  procedure  were  followed  here,  the  results  indicated  that  an 
average  step  size  somewhere  between  0.  07  and  0.  08  -  J  could  be  used. 

This  improvement  alone,  would  cut  the  calculation  time  by  a  factor  of 
approximately  two.  (A  short  test  case  indicated  this  to  be  feasible). 


Another  way  the  computation  time  could  be  decreased  would  be 
to  use  fewer  streamlines  throughout  the  flow  field.  In  the  present  computation, 
streamlines  were  added  at  the  shock  on  an  average  of  once  every  10  steps, 
so  that  a  total  of  21  points  were  specified  along  each  ray  in  the  final  plane. 

This  is  probably  more  than  enough  to  adequately  describe  the  flow  field  for 
most  ideal  gas  problems.  By  dropping  selected  streamlines  (e.  g.  ,  every 
other  one),  computing  time  would  be  reduced  because  one  needs  to  calculate 
fewer  points,  and  in  addition  the  axial  step  sizes  can  be  increased  as  a  direct 
consequence  of  the  larger  spacings  between  data  points.  Between  10  and  20 
points  on  each  ray  should  be  adequate  to  solve  steady  flow  problems  for  an 
ideal  gas. 


The  shock  wave  shape  and  some  typical  streamlines  are  shown  in 
Figure  20.  Since  the  scales  used  along  the  abscissa  and  ordinate  are  different, 
the  shapes  of  these  curves  are  somewhat  distorted;  however,  the  body  stream¬ 
line  does  correspond  to  a  15°  cone.  As  can  be  seen,  the  shock  shape  agrees 

very  well  with  that  obtained  from  the  two-independent  variable  characteristics 
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solution,  the  maximum  variation  amounting  to  about  0.  7%  at  x  =  7.  5.  The 
body  shape  and  streamline  patterns  are  also  as  would  be  expected.  Note  the 
manner  in  which  the  program  adds  streamlines  along  the  shock,  as  the 
calculations  proceed  in  the  axial  direction.  The  addition  of  these  extra 
streamlines  is  seen  to  be  absolutely  necessary  if  sufficient  accuracy  is  to 
be  maintained  throughout  the  flow  field. 

Results  for  the  flow  quantities  at  the  shock  and  body  in  one 
meridional  plane  are  presented  in  Figures  21  through  24.  In  all  cases,  the 
data  remained  axially-symmetric  to  at  least  four  significant  figures  over  the 
entire  range  of  the  computation.  Figure  21  shows  a  comparison  between  the 
pressures  obtained  from  the  present  program,  and  the  two  variable  program. 
The  results  are  in  good  agreement  and  indicate  that  the  truncation  errors 
are  approximately  the  same  for  both  methods  (order  of  the  step  size  squared). 
The  pressure  along  the  body  reaches  a  minimum  at  approximately  x  =  4.  25 
and  then  recompresses  toward  a  constant  sharp  cone  value  as  is  expected. 

The  total  variation  in  pressure  across  the  shock  layer  is  nearly  8%  greater 
using  the  three  variable  characteristics  scheme  than  with  the  two  variable 
method  at  an  axial  location  X  -  1.5.  However,  this  discrepancy  seems  to 
be  decreasing  and  as  cone  flow  is  reached,  it  should  nearly  disappear.  The 
velocity  /  2-D  flow  angle  Q  ,  and  temperature  T  are  shown  in  Figures  22, 

23  and  24  respectively,  compared  to  the  expected  axisymmetric  results. 

These  comparisons  again  clearly  show  the  accuracy  of  the  technique  with  the 
maximum  variations  amounting  to  less  than  1.  5%. 

Note  in  Figures  21  and  22  the  appearance  of  a  rather  severe 
expansion  for  the  3-D  solution  on  the  body  at  the  point  of  tangency  between 
the  spherical  nose  and  the  conical  afterbody  (  ,r=-0.  24)  compared  to  the 
axisymmetric  value.  The  over-expansion  is  approximately  10%  of  the 
expected  value,  and  is  propagated  throughout  the  flow  field  in  the  form  of 
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Figure  21  PRESSURE  DISTRIBUTION  ON  BLUNT  15°  CONE,  CC  =  0* 
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an  expansion  wave.  Figures  25,  26  and  27  show  the  propagation  of  this  wave 
along  rays  normal  to  the  body  axis  for  three  different  axial  locations.  As 
shown  in  Figure  j0  the  expansion  has  just  begun  to  form  ni  ir  the  body  at  an 
axial  location  of  x  =  -0.  231  and  the  pressure  throughout  most  of  the  shock 
layer  corresponds  to  the  expected  axisymmetric  value. 

4.  3.  2  Blunted  Cone  at  10°  Angle-of-Attack 

The  flow  field  about  the  same  15°  blunted  cone  was  computed  at 
10°  angle-of-attack  using  initially  16  meridional  rays  located  11.  33°  apart. 
Because  of  computing  time  limitations,  the  computations  were  halted  after 
an  X  coordinate  equal  to  3.  1  nose  radii  was  reached.  At  the  time  the 
computation  was  stopped,  it  was  proceeding  smoothly  and  could  have  been 
continued,  if  desired.  The  calculations  were  proceeding  at  the  rate  of 
approximately  one  nose  radius  per  hour  on  the  IBM  7094.  However,  the 
same  remarks  concerning  step  size  and  mesh  spacing  that  were  given  above 
for  the  zero  angle-of-attack  case,  would  also  apply  here.  As  an  example, 
the  average  step  size,  A  x  ,  used  in  the  calculation,  was  about  Ax  =0.03 
(nondimensionalized  with  respect  to  nose  radius)  while  a  short  run  indicated 
that  step  sizes  as  high  as  A%  =  0.  07  could  be  used  without  causing  instability 
for  at  least  20  axial  steps.  It  appears  that  optimizing  the  step  size,  perhaps 
using  the  method  discussed  previously,  is  required  in  order  to  improve  the 
efficiency  of  the  calculation. 

Another  way  the  running  times  could  be  reduced  would  be  to  use 
the  body  axis  as  the  X  -  axis  for  the  computation,  rather  than  a  wind - 
oriented  axis  (See  Figure  19)  as  was  utilized  here.  The  shift  to  the  new 
axis  could  be  accomplished  by  rotating  the  x-axis  to  coincide  with  the  body 
axis.  This  would  have  the  effect  of  decreasing  the  maximum  flow  angles 
and  as  a  result  increase  the  effective  step  size. 

The  technique  of  dropping  selected  "rings"  of  data  points  was 
successfully  tried  in  this  calculation.  When  10  extra  rings  had  been  added, 
for  a  total  of  19,  every  other  ring  throughout  the  field  was  dropped  from 
the  computation.  The  step  size  increased  by  a  factor  of  2  as  would  be 
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gure  26  PRESSURE  DISTRIBUTION  THROUGH  SHOCK  LAYER. BLUNT  I 
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Figure  28  FLOW  ANGLES  ALONG  SELECTED  STREAMLINES  FOR 
BLUNT  15°  CONE,  a=  0° 
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expected  and  there  was  no  noticeable  t  ifect  on  the  flow  properties  along 
streamlines.  The  computation  proceeded  for  0.  5  nose  radii,  during  which 
time,  another  3  rings  of  data  points  had  been  added  to  the  field.  The  overall 
computation  time  had  thus  been  decreased  by  a  factor  of  approximately  4  as 
compared  to  the  previous  rate.  The  results  of  the  calculation  are  shown  in 
Figures  29  through  36.  Figure  29  shows  a  side  view  of  the  body  and  the  shape 
of  the  shock  wave  on  the  most  windward  and  most  leeward  sides  of  the  cone. 
Also  shown  are  isobaric  curves  on  the  body  surface.  Note  the  region  of 
transition  from  axially  symmetric  flow  to  fully  three-dimensional  flow  on 
the  body  sui  ace.  The  pressures  on  the  windward  and  leeward  sides  of  the 
cone  differ  by  a  factor  of  nearly  five  at  the  last  axial  location  calculated. 

The  three-dimensionality  of  the  ttow  along  the  body  surface  can 
also  be  seen  in  Figure  30,  where  the  azimuthal  variation  in  pressure  around 
the  body  for  several  different  axial  stations  is  presented.  At  an  axial  location 
X  =  -0.  2,  the  pressure  is  still  constant  over  one  half  of  the  body  surface. 

By  the  time  station  X  =  0.  0  is  reached,  the  flow  is  completely  asymmetric, 
with  the  pressure  on  the  leeward  surface  having  fallen  by  a  factor  of  two. 

The  pressure  then  continues  to  decrease  around  the  cone  until  section 
X  -  1.0  is  reached.  Between  axial  stations  X  -  1 . 0  and  X=  2.0,  the 
pressure  near  the  windward  side  of  the  cone  has  increased  while  the  pressure 
over  the  remaining  portion  of  the  body  surface  continues  to  decrease.  When 
the  computation  was  stopped,  the  region  of  increasing  pressure  encompassed 
about  one-third  of  the  cone  surface,  and  was  progressing  toward  the  leeward 
side  as  expected.  The  cone  pressure  at  which  the  pressures  on  the  windward 
and  leeward  sides  should  stabilize  are  also  shown  and  appear  to  be  the  correct 
magnitude. 

The  axial  pressure  distribution  on  the  most  windward  generatrix 
of  the  cone  is  shown  in  Figure  31.  The  typical  overexpansion  is  evident  with 
the  minimum  occurring  in  this  case  at  approximately  one  nose  radius  from 
the  center  of  the  sphere.  As  shown,  the  pressure  is  approaching  the  equiva¬ 
lent  cone  pressure  and  should  be  constant  after  approximately  5  radii. 

The  distribution  of  pressure  across  the  shock  layer  is  shown  in 
Figures  32  and  35  in  planes  representing  four  different  axial  stations.  Note 
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that  the  body  shape  in  each  successive  plane  does  not  remain  circular. 

Again,  this  is  because  the  A-axis  is  parallel  to  the  free  stream  velocity 

vector.  The  results  presented  display  basically  the  same  pattern  as  those 

56 

which  occur  in  cross -sections  normal  to  the  body  axis.  After  the  non- 
axisymmetry  of  the  flow  has  been  established,  the  pressure  is  nearly  constant 
along  rays  approximately  normal  to  the  body  through  most  of  the  shock  layer. 
The  maximum  pressure  variation  through  the  shock  layer  is  approximately 
half  the  total  azimuthal  variation  on  the  body  at  X  =  3.124  with  the  largest 
gradients  existing  near  the  shock.  Note  also  that  the  shock  wave  remains 
nearly  circular  with  respect  to  the  wind  oriented  axis,  for  the  region  of  the 
flow  field  calculated.  The  scale  used  in  Figure  36  is  different  from  that 
used  in  the  three  previous  figures.  Hence,  although  the  body  diameter  and 
shock  layer  thickness  appear  to  be  smaller  than  those  at  station  X  -  2.  011 
they  are  actually  much  larger. 

The  pressure  distribution  obtained  along  the  body  surface  was 
smooth,  although  small  oscillations  similar  to  those  discussed  above  for  0° 
angle-of-attack  did  appear  along  rays  normal  to  the  body.  The  amplitudes 
of  the  oscillations  were  in  this  case  much  smaller  than  before,  and  disappeared 
when  the  expansion  wave  reached  the  shock.  The  reason  the  amplitudes  were 
smaller  can  probably  be  attributed  to  the  fact  that  smaller  axial  step  sizes 
were  taken  in  the  region  of  the  shoulder  between  the  nose  sphere  and  the  cone. 
Hence,  the  speculation  that  the  oscillations  can  be  completely  eliminated  if 
sufficiently  small  step  sizes  are  taken  in  this  region  appears  to  be  justified. 

The  pressure  distributions  shown  in  Figures  32  through  35  represent  the 
mean  values  of  the  oscillations.  A  typical  oscillation  around  the  mean  is 
show'n  in  Figure  34. 

The  intersections  of  various  streamlines  with  the  plane  X-  3.  124 
are  shown  in  Figure  36.  The  dashed  lines  represent  the  locations  of  the 
original  rays  from  the  body  to  the  shock,  as  they  would  appear  if  the  body 
and  flow  field  were  axisymmetric  with  respect  to  the  origin.  For  comparison, 
the  shap^  of  the  body  cross-section  as  it  appeared  in  the  initial  plane  is  also 
presented.  Each  curve  represents  the  locus  of  streamlines  which  were 
initially  located  on  the  same  ray.  The  streamlines  near  the  body  are  "spinning" 
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Figure  34  PRESSURE  DISTRIBUTION  ACROSS  SHOCK  LAYER  FOR  SPHERICALLY  BLUNTED 
CONE,  p  =  p'//2.'U£ 
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Figure  35  PRESSURE  DISTRIBUTION  ACROSS  SHOCK  LAYER  FOR  SPHERICALLY  BLUNTED 


toward  the  leeward  side,  with  the  maximum  displacements  away  from  the 
initial  ray  having  occurred  on  the  body  in  the  region  where  the  pressure 
gradients  are  greatest,  as  would  be  expected.  Near  the  shock  however,  the 
displacements  are  slight,  indicating  again  that  the  flow  in  this  region  has 
remained  nearly  axisymmetric.  The  fact  that  the  streamlines  flow  toward 
the  leeward  side  indicates  that  it  will  be  necessary  to  include  a  subroutine 
in  the  program  which  adds  new  streamlines  through  the  shock  layer  on  the 
windward  side  as  they  are  needed  for  sufficient  resolution  of  the  flow  field. 

One  convenient  way  to  allow  for  this,  would  be  to  choose  more  initial  rays 
on  the  windward  side  of  the  body,  and  use  unequal  meridional  spacings  for 
the  rays  in  the  initial  plane.  The  program  is  capable  of  handling  this  situation. 
Also,  it  will  be  necessary  to  drop  streamlines  as  they  crowd  together  on  the 
leeward  side.  Probably  the  best  way  to  accomplish  all  of  these  objectives 
would  be  to  realign  the  mesh  along  body  normals  by  interpolation  at  selected 
axial  stations.  Although  this  procedure  involves  shifting  from  one  set  of 
streamlines  to  another  midway  through  the  calculation  scheme,  it  can  be 
accomplished  with  little  extra  programming  effort,  since  the  interpolation 
routine  is  already  available  in  the  program.  The  percentage  increase  in  the 
overall  computation  time  should  be  small,  since  the  mesh  points  need  be 
realigned  in  this  manner  perhaps  only  three  or  four  times  throughout  a 
particular  flow  field  calculation. 


4.  3.  3 


Elliptical  Afterbod\ 


The  third  example  consists  of  the  flow  field  around  a  spherically 
blunted  elliptical  alterbody  at  zero  angle-of-attack.  The  locus  of  the  semi¬ 
major  axes  of  the  ellipses  in  cross-sections  normal  to  the  body  axis  is  a 
straight  line  inclined  at  15°  to  the  body  axis,  while  the  locus  of  the  semi¬ 
minor  axes  is  a  straight  line  inclined  at  5°  to  the  axis.  Thus,  cross-sections 
of  the  body  are  elliptical  with  eccentricities  varying  as  a  function  ef  the 
axial  coordinate  A  .  The  equation  of  the  body  is  written  in  three  separate 
sections.  First,  the  spherical  nose; 
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Second,  the  transition  region  through  which  the  locus  of  semi-major  axes 
is  a  straight  line,  while  the  locus  of  semi-minor  axes  is  a  circle; 
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where  a(x )  -  ( 2.  X  -  X  *)  ,7*h  t  !8  <  x*  .11*9+  (130) 

b  (*)  =  •  7d  73  +  .  ZC,  8  X  (131) 


and  finally,  the  afterbody  region  where  the  loci  of  the  semi-major  and  minor 
axes  are  both  straight  lines  and 
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where , 

Cl(x)  =  .  9!  Z#+ +.027+9  X  (133) 


h>  (x)  =  .  2  +  8*  +  .  76  73  (134) 
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The  resv.its  of  the  computation  are  given  in  Figures  37  through 
40.  In  this  case,  a  total  of  9  meridional  .  s  containing  9  points  each 
was  used  in  the  initial  plane.  The  ray  spacings  were  constant  and  equal 
to  11.  25* .  It  required  nearly  90  minuter  on  an  IBM  7094  to  calculate  3 
nose  radii  using  an  average  step  of  0.03  radii.  Seven  rings  of  data  points 
were  added  at  the  shock  throughout  the  entire  computation,  while  no  point 
dropping  was  employed.  Again,  since  rather  large  step  sizes  (0.08)  were 
used  in  the  transition  region  near  the  nose,  oscillations  similar  to  those 
discussed  above  did  develop  in  the  field.  The  curves  shown  represent 
the  mean  values  of  the  oscillations.  The  amplitudes  of  these  oscillations 
were  again  dec reas ing  as  the  compu  tation  ;  ,  occeded  downstream.  However, 
since  the  initial  expansion  wave  bad  not  yet  reached  the  shock,  the  oscillations 
were  still  apparent.  A  typical  oscillation  around  an  isobar  in  the  final  data 
plane  calculated  ip  shown  in  Figure  40.  In  this  case,  the  expansion  wave  is 
just  c-ver  one-half  the  distance  through  the  shock  layer. 
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The  pressure  distribution  along  the  body  is  shown  in  Figure  37 
in  the  form  of  isobaric  curves.  The  locus  of  the  semi-minor  axes  coincides 
with  the  X-axis,  while  the  locus  of  the  semi-major  axes  is  shown  in  the 
figure.  The  asymmetry  of  the  distribution  is  apparent,  with  a  pressure 
decrease  of  30%  occurring  around  the  body  in  the  last  plane  calculated 
(X  =  3.  0).  The  "humps"  which  occur  along  the  isobars  near  the  nose  are 
again  caused  by  the  initial  overexpansion  in  the  transition  region.  They 
could  be  eliminated  by  taking  smaller  step  sizes  in  this  region. 

The  pressure  distributions  through  the  shock  layer  are  presented 
in  Figures  38,  39  and  40  for  cross-sections  at  X  =  1 . 07,  2.  01  and  3.  02, 
respectively.  The  shock  wave  and  body  shapes  are  also  shown  in  each  plane. 
Note  that  the  shock  is  still  axisymmetric  in  the  last  plane  calculated.  Also, 
the  bluntness  factor  of  the  body  ellipse  in  the  last  plane  calculated  is  only 
1.  2  compared  with  the  asymptotic  value  of  2.  2  for  this  particular  body.  Thus, 
the  problem  is  not  nearly  as  three-dimensional  as  was  intially  anticipated. 
Computer  time  limitations  prevented  carrying  the  solution  out  farther  in 
the  axial  direction. 

In  the  last  plane,  nearly  all  of  the  shock  layer  has  an  asymmetric 
pressure  distribution,  whereas,  at  X  ~  1.07,  the  distribution  is  still  axi  - 
symmetric  through  almost  one- half  of  the  layer.  Along  the  low  pressure 
edge  of  the  body,  the  shock  layer  thickness  is  1  .  3  as  compared  with  the 
value  of  1 . 0  along  the  high  pressure  edge  in  this  final  plane.  The  shock 
layer  thickness  is  thus  over  three  times  as  large  as  the  initial  value  of 
0.  37.  Generally,  the  pressure  distributions  correspond  to  what  would  be 
expected  for  the  elliptical  afterbody,  and  except  for  the  oscillations  which 
did  arise  through  the  field,  the  results  of  the  run  seem  to  again  verify  the 
fact  that  the  present  version  of  the  program  is  working  correctly. 


Figure  39  PRESSURE  DISTRIBUTION  ACROSv  SHOCK  LAYER  FOR  SPHERICALLY  BLUNTED 
ELLIPTICAL  AFTERBODY,  <*  =  O*,  X  =  2.014 
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Figure  40  PRESSURE  DISTRIBUTION  ACROSS  SHOCK  LAYER  FOR  SPHERICALLY  BLUNTED 
ELLIPTICAL  AFTERBODY,  Ct  =  0°,  X  =  3.018 
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Section  5 


CONCLUSIONS 

Many  three-dimensional  steady  flow  fields  can  be  calculated 
utilizing  the  numerical  procedures  presented  in  this  report;  however,  the 
computer  times  are  quite  long.  Extension  of  the  program  developed  here'  to 
include  physico-chemical  models  more  complex  than  the  perfect  gas  can  also 
be  readily  accomplished.  As  a  matter  of  fact,  the  scheme  developed  is 
especially  applicable  to  solving  chemically  frozen  and  finite-rate  reacting 
flows.  However,  the  economic  practicality  of  such  calculations  is  considered 
to  be  marginal  with  presently  available  equipment. 

On  the  next  generation  of  computers  these  problems  should  be 
feasible.  The  three-dimensional  steady  flow  program  developed  here  involve 
large  amounts  of  input  and  output  data.  Although  the  abundance  of  output 
data  can  be  controlled  by  printing  out  only  a  few  selected  planes  of  data  during 
a  computation,  the  tasks  of  handling  and  presenting  even  these  data  are  still 
tedious.  The  development  of  internal  computer  methods  for  quickly  handling 
and  displaying  the  output  data  graphically  or  in  some  other  convenient  form, 
is  thus  necessary. 

Several  basic  conclusions  can  be  stated  regarding  the  general 
numerical  scheme.  First,  the  Courant-Friedrichs-Lewy  (C-F-L)  stability 
criterion  proved  to  be  a  sufficient  condition  to  insure  stability  of  the  integra¬ 
tion  procedure.  It  should  be  noted,  however,  that  with  the  fitting  procedures 
utilized  here,  strict  adherence  to  the  C-F-L  condition  is  probably  not  always 
required.  Hence,  it  appears  possible  to  enlarge  the  domain  of  dependence  of 
the  difference  equations  slightly  by  extrapolation  using  the  polynomial  surface 
fits  without  affecting  the  accuracy  of  the  solution.  The  total  computing  time 
per  problem  would  in  turn  be  decreased.  Second,  for  accuracy  and  stability, 
higher  order  fitting  procedures  than  linear  are  necessary.  The  quadratic  fits 
used  here  are  expected  to  be  sufficiently  accurate  for  most  flow  problems. 
Finally,  in  order  to  reduce  truncation  errors  to  second  order  in  the  step  size, 
it  is  necessary  to  add  extra  bicharacteristics  in  addition  to  those  that  would 
be  required  for  a  simplical  mesh,  The  results  obtained  from  the  program 
utilizing  the  extra  bicharacteristics  compared  closely  with  those  obtained 
from  a  two  independent  variable  method  of  characteristics  solution. 
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The  greatest  advantage  to  be  derived  from  the  procedure 
developed  here  lies  in  the  fact  that  for  the  most  part  (except  when  mesh 
"squaring"  is  used)  the  integrations  proceed  along  the  same  streamlines. 

This  greatly  simplifies  inhomogeneous -frozen  and  nonequilibrium  flow  field 
calculations  by  significantly  reducing  the  number  of  interpolations.  Having 
the  output  data  available  along  streamlines  is  also  useful  for  visualizing 
the  flow  field  associated  with  a  particular  body,  and  for  applying  the  pressure 
distribution  obtained  from  the  program  to  calculate  three-dimensional  boundary 
layers . 

For  the  Cauchy  problem  considered,  a  large  amount  of  initial 
data  is  required  to  calculate  a  given  flow  field.  For  a  general  three-dimen¬ 
sional  body,  this  data  must  be  supplied  at  points  on  an  initial  plane  (or  surface) 
located  entirely  in  the  supersonic  portion  of  the  flow  field.  Thus,  what  is  in 
general  required  is  some  method  for  obtaining  asymmetric  inputs  on  the 
initial  plane.  Unfortunately,  such  a  method  is  not  yet  available,  and,  for 
the  problems  considered,  it  is  necessary  to  limit  the  geometries  to  spherically 
capped  afterbodies,  for  which  axisymmetric  inputs  can  be  utilized.  When 
methods  for  calculating  non-axisymmetric  blunt  body  flows  become 
available,  they  can  be  used  to  supply  inputs  to  the  steady  flow  method  of 
characteristics  scheme  utilized  here.  One  possible  way  to  obtain  asymmetric 
inputs  would  be  to  generalize  the  present  program  to  include  unsteady  flow. 
Starting  from  one  known  steady  blunt  body  solution,  another  steady  flow  solu¬ 
tion  could  be  obtained  asymptotically  in  time,  after  the  nose  cap  has  undergone 
certain  prescribed  motions  (or  shape  changes).  Thus,  using  a  combination 
of  three-dimensional  steady  and  unsteady  program,  one  could  analyze  rather 
general  body  shapes  and  motions. 
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Appendix  A 

PROGRAM  DETAILS  AND  OPERATING  PROCEDURES 


This  appendix  presents  the  information  necessary  for  running 
the  program  and  interputing  the  output.  Paragraph  A.  1  contains  information 
on  the  program  limits  and  its  application  and  tape  requirements.  Paragraph 
A.  2  gives  detailed  description  of  the  input  data  and  its  preparation  for 
starting  or  restarting  a  given  computational  problem.  Also  included  here 
are  prepared  input  format  sheets  for  ease  in  punching  the  input  data  cards. 

The  output  from  the  program  is  discussed  in  Paragraph  A.  3.  Details  con¬ 
cerning  the  data  storage  are  given  in  Paragraph  A.  4.  A  description  of 
each  subroutine  is  included  in  Paragraph  A.  5.  A  correspondence  table 
listing  the  correspondence  between  the  FORTRAN  and  the  report  symbols 
is  included  in  Paragraph  A.  6.  Finally,  Paragraph  A.  7  includes  logical 
flow  charts  describing  the  program. 

A.  1  DESCRIPTION  OF  DIGITAL  COMPUTER  PROGRAM 

The  computer  program  described  here  is  an  extension  and 
modification  of  that  developed  in  Reference  59.  It  was  originally  intended 
that  the  final  program  would  be  capable  of  computing  supersonic  flows 
arour J  fairly  general  bodies  including  exactly  the  effects  of  chemical 
reactions  and  vibrational  relaxations.  However,  experience  with  the 
ideal  gas  three  variable  program  and  axisymmetric  and  two-dimensional 
solutions  (see  Reference  60)  indicated  that  such  computations  were  beyond 
the  capabilities  of  present  computer  systems  \chen  reasonable  running 
times  and  machine  storage  requirements  were  taken  into  account.  Therefore, 
in  the  interest  of  conserving  machine  time  while  at  the  same  time  obtaining 
a  reasonably  close  approximation  to  the  flow  field,  it  was  decided  that  the 
next  logical  and  least  costly  extension  to  the  ideal  gas  program  would  be 
to  add  the  frozen-inhomogeneous  thermo-chemical  model  described  above. 
This  modification  was  introduced  into  the  present  version  of  the  program. 

The  pressure  distribution  thus  obtained  corresponds  to  the  frozen  solution. 


1  15 


0 


It  may  then  be  utilized  to  obtain  a  more  complete  description  of  the  flow 
field  by  using  a  one  dimensional  stream  tube  integration  program  aiong 
each  of  the  streamlines  followed  in  the  characteristics  procedure.  To 
facilitate  the  computations  for  the  nonequilibrium  solution,  the  streamline 
locations  and  pressures  for  each  step  are  printed  out  on  binary  tape  so  that 
they  may  then  be  easily  used  by  a  suitable  stream  tube  integration  program. 
Actually  such  a  computation  has  been  performed  but  the  program  required  is 
not  yet  completed  and  is  not  being  reported  here.  Since  the  present  version 
of  the  program  is  considerably  different  than  the  original  ideal  gas  version, 
the  solutions  described  above  and  in  Reference  59  were  recomputed  for  an 
ideal  gas.  In  this  manner  many  program  errors  were  uncovered  and 
corrected.  The  program  was  also  run  using  the  frozen-inhomogeneous 
option  for  a  6  inch  nose  radius  spherically  capped  cone  at  zero  angle  of 
attack.  The  results  indicated  that  this  option  is  working  correctly. 
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A.  1. 1 


General  Discussion 


The  greatest  part  of  the  computations  involved  in  determining 
the  flow  properties  in  the  shock  layer  about  a  three-dimensional  body 
are  contained  in  the  subroutines  which  calculate  body  points,  field  points 
and  shock  points  in  the  manner  set  forth  in  Section  3.  All  other  sub¬ 
routines  (incl  rng  the  main  program)  are  merely  subordinate  routines 
to  these  three  and  a.e  used  to  shift  required  data  around  in  storage, 
perform  supporting  calculations  and  handle  input/output  procedures. 

The  main  program  controls  the  overall  logic  of  the  scheme  by 
first  choosing  an  initial  point  located  in  the  field  or  on  the  body  or  shock 
surfaces,  calling  the  respective  principal  subroutine  and  accepting  the 
results  for  the  newly  calculated  point.  It  then  stores  these  results  as  a 
point  the  new  data  plane  and  proceeds  on  to  another  initia*  point.  When 
all  the  points  on  the  initial  data  plane  are  thus  determined  the  newly 
calculated  ciata  plane  then  serves  as  the  next  initial  plane,  and  the 
calculations  are  continued.  New  data  planes  normal  to  the  X -direction 
are  continually  calculated  until  either  a  predetermined  coordinate  is 
reached,  or  a  given  maximum  number  of  planes  are  so  determined  --  at 
which  time  the  computations  are  halted.  The  relationships  among  the 
various  subroutines  used  and  the  main  program  are  shown  in  terms  of  a 
general  calling  sequence  chart  in  Figure  A-l.  A  detailed  description  of 
each  of  the  subroutines  will  be  given  in  Paragraph  A.  3. 

A.  1.2  Program  Versions 

Due  to  the  length  of  the  program,  it  is  necessary  to  employ 
either  the  CHAIN  feature  of  the  7044  operating  system  or  an  overlay 
structure  on  the  7094  (see  Figure  A-2).  The  differences  between  the  two 
operating  systems  has  necessitated  the  preparation  of  two  closely  similar 
versions  of  the  program.  The  version  described  here  is  that  suitable  for 
the  7094.  For  the  7044  system  the  subroutine  INPUT  is  treated  as  a  separate 
link  in  the  CHAIN, 


S 
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Figure  A- 2  OVERLAY  STRUCTURE  FOR  IBM  7094 


A.  1.3  Tape  Requirements 

(a)  TAPE  2  binary  tape.  Data  from  an  integration 

plane  are  stored  on  tape  2  for  every  N 
planes  computed  where  N  is  an  input  constant. 
This  tape  is  used  for  restarting  a  computation. 

(b)  TAPE  3  binary  tape.  Common  data  are  stored  on 

tape  3.  This  may  be  used  to  provide  inputs 
to  a  finite  rate  stream  tube  integration  program. 

(c)  TAPE  4  binary  tape.  Streamline  data  which  can  be 

used  for  reacting  stream  tube  integration 
scheme  is  stored  on  binary  tape  4. 

(d)  TAPE  5  System  input  tape. 

(e)  TAPE  6  System  output  tape. 
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DESCRIPTION  AND  PREPARATION  OF  INPUT  DATA 


A.  2 


As  was  discussed  in  Paragraph  3.  5,  a  plane  (or  surface)  of  inputs 
is  required  to  any  three-dimensional  characteristics  program.  Obtaining 
data  for  such  an  initial  plane  is  sometimes  difficult  owing  to  the  lack  of 
suitable  three-dimensional  blunt-bodv  solutions.  However,  for  many  config¬ 
urations  having  spherical  noses,  it  is  possible  to  assume  axisymmetric  flow 
up  to  some  axial  station.  For  such  cases,  the  program  described  in 
Reference  oO  may  be  modified  to  punch  out  most  of  the  input  cards  needed 
for  the  present  program.  The  necessary  changes  are  detailed  in  Appendix  C. 


In  this  report,  the  inputs  were  obtained  by  combining  an  inverse 

39  35 

technique  with  an  axisymmetric  method  of  characteristics  solution  .  The 

inverse  scheme  of  Garr  and  Marrone  was  used  to  obtain  several  data  points 
between  the  shock  and  the  body  on  a  line  normal  to  the  shock  and  located 
slightly  downstream  of  the  sonic  line.  The  two  variable  method  of  character¬ 
istics  program  was  then  utilized  to  extend  this  solution  a  short  distance  down¬ 
stream,  until  a  sufficient  number  of  initial  data  points  could  be  determined 
on  a  line  normal  to  the  free  stream  wind  direction  (see  Figure  19).  A  plane 
normal  to  the  wind  direction  containing  this  line  is  chosen  to  be  the  initial 
data  plane  for  the  three-dimensional  computation.  Since  the  body  is  initially 
axisymmetric  (e.  g.  ,  a  spherical  nose  cap)  it  suffices  to  obtain  the  input  along 
one  line  only.  Flow  properties  at  initial  data  points  throughout  the  initial 
plane  are  obtained  by  rotating  the  initial  data  line  around  the  x  axis  (e.  g.  , 
the  wind  axis)  in  chosen  increments  in  <f>  ,  the  angle  between  any  azimuthal 

plane  and  the  x  y  plane.  At  each  point  chosen  on  the  initial  data  line,  the 
radius  f  ,  the  pressure  P  ,  the  velocity  ^  .  and  the  two-dimensional  flow 
angle  X  are  required  as  inputs  to  the  program.  In  addition,  the  two- 
dimensional  shock  angle  a~  ,  the  angle-of-attack  cc.  ,  the  gas  constant  R  , 
and  various  free  stream  and  nondimensionali zing  parameters  are  needed. 

All  input  parameters  are  dimensioned  such  that  the  following  form  of  the 
equation  of  state  holds 


F  9 1° 


/  £' 

AM" 


T' 


(A-  1 ) 
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for  whatever  units  are  chosen.  Energies  are  made  nondimensional  with 
relations  of  the  form 


(A-2) 


while  for  pressure 


(A-  3) 


is  used.  Here,  the  subscript  o®  refers  to  free  stream  conditions  while 
the  primes  denote  dimensional  quantities.  In  the  computer  program  the 
number  0  (zero)  is  used  to  designate  free  stream  conditions  and  is  inter¬ 
changeable  with  oo  .  As  long  as  Equation  (A-l)  through  (A- 3)  are  satisfied 
any  system  of  units  may  be  employed.  A  complete  listing  of  the  required 
inputs  and  associated  formats  is  given  in  Paragraph  A.  5. 

It  should  again  be  mentioned  here  that  this  method  of  starting 
the  three-dimensional  characteristics  solution  from  a  plane  of  axisymmetric 
inputs,  is  valid  only  when  the  geometry  downstream  of  the  sphere  is  such 
that  the  sonic  point  is  located  on  the  sphere,  for  the  geometries  and  angles 
of  attack  being  considered. 


A.  2.  1  Description  of  Input  Data 

Input  to  the  program  is  provided  on  data  cards  if  the  problem  to  be 
computed  is  a  new  case.  If  the  problem  is  a  restart  of  a  previous  incomplete 
run,  input  data  are  provided  on  cards  and  on  a  binary  tape.  A  description  of 
the  format  and  order  of  the  input  data  cards  follows  is  shown  in  Table  A-l. 
The  input  formats  have  been  simplified  as  much  as  possible,  with  1615 
being  used  for  fixed  point  data  and  5E1  5.  7  being  used  for  floating  point  as 
they  are  consistent  with  the  formulae  in  Paragraph  A.  2.  The  steps  of  this 
input  procedure  are  keyed  to  the  blanks  in  the  accompanying  INPUT  FORMAT 
(See  Figure  A-3).  The  sign  of  the  exponent  is  located  by  a  small  s.  Enter 
decimal  points  and  right-adjust  integer  numbers. 
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TABLE  A-I 


TABLE  OF  INPUT  DATA  CARDS 


Ca  rd 
Number 

Column 

Number 

MNEMONIC 

Item  Entered 

Number 

R  epr  esentation 

1 

1  -  5 

NR  UN 

Control  labeling 
run  number. 

Any  integer 
number  may  be 
used. 

Integer 

15 

2** 

1  -  5 

ITAPE 

Control  constant  - 
If  ITAPE  is  not 
equal  to  zero, 
data  store  is  on 
binary  tape. 

Integer 

15 

6  -  10 

N  P  LIN 

Control  constant  - 
When  ITAPE  f  0, 
NPLIN  defines  the 
number  of  the  data 
plane  from  which 
the  computations 
are  to  be  re¬ 
started. 

Integer 

15 

11  -  15 

IREC4N 

Control  constant  - 
defines  number  of 
records  to  be 
skipped  on  tape  4 
when  ITAPE  £  0. 

Integer 

15 

3 

1  -  5 

NSPEC 

Number  of  species 
included  in  gas 
model  description 

2  i  NSPEC  £  8 

Integer 

15 

6  -  10 

NR 

Not  used  in  frozen 
program  -  Number 
of  reactions  in¬ 
cluded  in  non¬ 
equilibrium  option 

Integer 

15 

11-15 

NC 

Number  of  elements  Integer 

15 

'"'"Note:  If  ITAPE  £  0,  the  only  remaining  data  cards  are  3*through  6* 
below;  otherwise  cards  3  through  6*  are  used 
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Table  A-I  (Cont.  ) 

TABLE  OF  INPUT  DATA  CARDS 


Card 

Number 


Colum.i  Number 

Number  MNEMONIC  Item  Entered  Representation 


16-20  NVIB 


21-25  MODEL 


26  -  30  NP180 


31-35  IPT5 


36  -  40  NPTS 


41  -  45  NROTAT 


46  -  50  MDELAY 


Number  of  vibra¬ 
tors  not  in  vibra¬ 
tional  equilibrium 

0  ^  NVIB  -  2 

Integer 

15 

0  -  frozen  solution 
is  computed 

Integer 

15 

i  -  ideal  gas  model 
is  assumed 

C  -  90°  section 
of  flow  field  is 
computed 

Integer 

15 

1  -  180°  section 
is  computed 

Number  of  points 
on  initial  data 
input  ray 

3  &IPTS  £  21 

Integer 

15 

Number  of  rays  in 
initial  data  plane 

3  <  NPTS  s=  29 

Integer 

15 

Control  constant  - 
Used  when  initial 
coordinate  system 
is  to  be  rotated  to 
body  fixed  system 

Integer 

15 

Control  constant  - 
Defines  number  of 
rings  wiiich  are 
used  in  surface 
fitting  procedure. 

Set  MDELAY  =  3 
always . 

Integer 

15 
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Table  A-I  (Cont.  ) 

TABLE  OF  INPUT  DATA  CARDS 


Ca  rd 
N  umber 


Column 

Number 


MNEMONIC  Item  Entered  Re 


4 


5 


6 


1  -  5  for  J  =  1  NS  PE C ( J ) 

6  -  1 0  f  o  r  J  =  2 
11  -  1  5  for  J  •-  3 
16  -  20  for  J  =  4 
21  -  25  for  J  =  5 
26  -  30  for  J  -  6 
31  -  35forJ  =  7 
36  -  40  for  J  =  8 


1  -  15 

16  -  30 

1  -  15 

16  -  30 

31  -  45 

46  -  60 

ol  -  75 


SIGMA 

SIGMA  1 

RH0 

uo 

R0 

TO 

W  MW  0 


Define  which  species 
are  used  in  model 
description.  Species 
are  tagged  a.s  follows 

NOS  PEC  =  1  for  O 
NOSPEC  -  2  for  N 
NOS  PEC  =  3  for  c- 
NOSPEC  =  4  for  Ar 
NOSPEC  =  5  for  Oz 

NOSPEC  -  6  for  N£ 
NOSPEC  =  7  for  NO 
NOSPEC  =  8  for  NO  + 

2-D  shock  angle  at 
shock  point  on  initial 
ray  (degrees) 

2-D  shock  angle  at 
a  shock  point  slightly 
downstream  of  initial 
ray 

F  ree  stream  density- 
(gms/cm*)  or 
(slugs  / ft5) 

F  ree  stream  velocity 
(cm/sec)  or  (ft/sec) 

Universal  gas 
constant  (ergs/mole  °K) 
or  (ft  lb/mole  °R) 

F ree  stream 
temperature 
( °  K)  or  ( °R) 

F  ree  stream  molecular 
weight 

D  /IC 

/  9"  o-  — 5  \ 
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Number 

presentation 

Integer 

15 


Floating 
El  5.  7 


Floating 
El  5.  7 


Floating 

E15.7 


F  loating 
El  5.  7 

Floating 
El  5.  7 


F  loating 
El  5.  7 


F  loating 
El  5.  7 


0 


Table  A-l  (Gont.  ) 

TABLE  OF  INPUT  DATA  CARDS 


Card 

Number 


1  -  15 

DIML 

Reference  length 
(usually  bow  shock 
radius)  (cm)  or  (ft) 

F  loating 
El  5.  7 

16  -  30 

HINF 

Total  enthalpy 
(nondimensional) 

F  loating 
El  5.  7 

31  -  45 

XO 

x-coordinate  of 
initial  plane  -  Note: 
x,  y,  z  system  is 
referenced  at  center 
of  spherical  nose 
(nondimensional) 

^  *0 

F  loating 
El  5.  7 

46  t*-  60 

THETFS 

F  ree  stream  flow 
angle  (degrees) 

F  loating 
El  5.  7 

61  -  75 

PSIFS 

F  ree  stream  flow 
angle  (degrees) 

(Note:  always  enter 

0. 0  here) 

Floating 
E 1  5.  7 

1  -  15 

BIGAM 

F  ree  stream  specific 
heat  ratio 

F  loating 
E15.7 

16  -  30 

EPSINF 

(1) 

F  ree  stream  vibra¬ 
tional  energy  for 
fi  rst  vibrator 
(Note:  If  NVIB  =  0, 
EPSINF  (1)  =  0.0) 
(nondimensional) 

/r'tL 

F  loating 
E 1  5.  7 

31  -  45 

EPSINF 

(2) 

F  ree  stream  vibra¬ 
tional  energy  for 
second  vibrator 
(Note:  If  NVIB  -  0 

Floating 
El  5.  7 

or  1 ,  EPSINF(2)  =  0.  0) 
(Nondimensional) 

6-  <x>x.  ~  /R'^~o 


Column  Number 

Number  MNEMONIC  Item  Entered _  Representation 
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TABLE  A-I  (Cont.  ) 
TABLE  OF  INPUT  DATA*  CARDS 


Ca  rd 
Number 


o 

/ 


10 


1  1 


12 


Column  Number 

Number  MNEMONIC  Item  Entered  Representation 


46  -  60  PGRAD  Estimated  pressure  Floating 

gradient  along  El  5.  7 

streamlines 
0.  95£  PGRAD^  1. 05 
(Expansion)  (Compression) 


61-75  AXIS 


1-15  ATTACK 


16-30  ERROR 


1  -  15  for  J  s  1  GAMAM(J) 
16  -  30  for  J  =  2 
31  -  45  for  J  =  $ 

46  -  60  for  J  =  4 
61  -  75  for  J  =  5 


1  -  1  5  for  J  =  6  GAMA.M(J) 
16  -  30  for  J  =  7 
31  -  45  for  J  =  8 


1-15  ALPHA  (I) 


16  -  30  R(I) 


Angle  final  x-axis 
makes  with  body 
axis  after  rotation 
(degrees) 

Angle -of -attack 
(degrees) 


Floating 

E15.7 


F  loating 
E 1  5.  7 


Error  test  stop  for  Floating 

internal  program  El  5.  7 

iterations  (usually 
error  =  1 . 0  x  10"^) 


Free  stream  concen-  Floating 

tration  for  species  El  5.  7 

(J)  in  order  specified 
by  NOSPEC(J)  (see 
card  4)  (Note:  If  J 
>  NSPEC,  GAMAM(J) 

=  0.0).  (Nondimensional) 

(  AiU^wo/es 


Free  stream  concen-  Floating 

trations  for  species  El  5.  7 

(J)  (6  -  J  ~  8)  (see 
details  under  card  10) 


2-D  flow  angle  (  ac  )  Floating 
for  input  c.ata  point  El  5.  7 

(I)  (radians) 


2-D  radius  to  point  Floating 

(I)  on  initial  line  El  5.  7 

(nondimensional) 

“  /Vr  /dim l 
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TABLE  A-I  (Cont.  ) 

TABLE  OF  INPUT  DATA  CARDS 


Card 

lumber 

Column 

Number 

MNEMONIC 

Item  Entered  R 

Number 

epresentation 

31  -  45 

P(I) 

Pressure  for  point 
(I)  on  initial  line, 
(nondimensional) 

Floating 

E  1 5.  7 

46  -60 

Q(I) 

Velocity  for  point 
(I)  on  initial  line 
(nondimensional) 
o*  =  r  /Vac. 

Floating 

El  5.  7 

61  -  75 

T  (I) 

Temperature  for 
point  (I)  on  initial 
line  (nondimensional) 

-  TsAL 

F  loating 

E  1  5.  7 

13 

1  -  15 

RHO  (I) 

Density  for  initial 
point  (nondimensional) 

F  loating 

E  1  5.  7 

16  -  30 

XMW  (I) 

Molecular  weight 
for  initial  point  (I) 
(nondimensional) 

MWj  *  MWj  /HU/X, 

Floating 

E 1  5.  7 

14 

1  -  15 

EPSI(I,  1) 

Vibrational  energy  of 
first  vibrator  for 
point  I. 

(nondimensional) 

€>>i  ’rio, 

(EPSI(I,  1)  =  0.  0  if 

NVIB  =  0) 

Floating 

E15.7 

16  -  30 

EPSI(I,  2) 

Vibrational  energy  of 
second  vibrator  for 

Floating 

F,  1 5.  7 

point  (I). 
(nondimensional) 

(EPSI(I,  2)  =  0.  0  if 
NVIB  -  1) 
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TABLE  A-I  (Cont.  ) 
TABLE  OF  INPUT  DATA  CARDS 


Card 

Column 

Number 

Number 

Number 

MNEMONIC 

Item  Entered  R 

epresentation 

15 

1 

-  15, 

J  =  1 

GAMMA  1(1,  J) 

Concentration  of 

Floating 

16 

-  30, 

J  =  2 

species  (J)  for  point 

El  5.  7 

31 

-  45, 

J  =  3 

(I)  on  initial  data 

46 

-  60, 

J  =  4 

line,  (nondimensional) 

61 

-  75, 

J  =  5 

f  _  y  '  . ..  /  ^o.«5  j 

JnTi 

(if  J  NS  PEC  =  0.  0) 

16 

1 

-  15, 

J  =  6 

GAMMAIfl,  J) 

See  description  on 

Floating 

16 

-  30, 

J  =  7 

card  15 

E15.7 

31  -  45,  J  =  8 

Note:  Cards  12  through  16  are  repeated  for  each  point  on  the  initial  line 

from  the  body  (I  =  1 )  to  the  shock  (I  =  IPTS)  until  a  total  of  1  1  +  5IPTS 
cards  have  been  added.  This  completes  the  specification  of  the  initial 
ray  input  data. 


12 

+ 

5IPTS 

1 

- 

15, 

N 

2 

OMEGA(N) 

Azimuthal  angle  of 

Floating 

16 

- 

30, 

N 

= 

3 

ray  (N)  in  initial 

E  1  5.  7 

31 

- 

45, 

N 

= 

4 

plane  (degrees) 

46 

- 

60, 

N 

= 

5 

(if  N  >  NMAX, 

61 

- 

75, 

N 

6 

OMEGA(N)  =  0.  0) 

13 

+ 

5IPTS 

1 

- 

15, 

N 

= 

7 

OMEGA(N) 

See  card 

Floating 

16 

- 

30, 

N 

8 

12  +  5IPTS 

El  5.  7 

31 

- 

45, 

N 

= 

9 

46 

- 

60, 

N 

r 

10 

61 

- 

75, 

N 

1  1 

14 

+ 

5IPTS 

1 

- 

15, 

N 

= 

12 

OMEGA  (N) 

See  card 

F  loating 

16 

- 

30, 

N 

= 

13 

1.1  +  5IPTS 

E  1  5.  7 

31 

- 

45, 

N 

r 

14 

46 

- 

60, 

N 

= 

1  5 

61 

- 

75, 

N 

- 

16 

15 

+ 

5IPTS 

l 

- 

15, 

N 

= 

17 

OMEGA(N) 

See  card 

F  loating 

16 

- 

30, 

N 

= 

18 

12  +  5IPTS 

El  5.7 

31 

- 

4  5, 

N 

19 

46 

- 

60, 

N 

= 

20 

61 

- 

75, 

N 

21 
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TABLE  A  - 1  (Cont.  ) 

TABLE  OF  INPUT  DATA  CARDS 


Ca  rd 
Number 

Column 
Numbe  r 

MNEMONIC 

Item  Entered 

Numbe  r 

Re  pres  entatiori 

If  +  5IPTS 

1 

-  15, 

N  -  22 

OMEGA(N) 

See  card 

F  loating 

16 

-  30, 

N  =  23 

12  +  5IPIS 

El  5.  7 

31 

-  45, 

N  -  24 

46 

-  60, 

N  =  25 

61 

-  75, 

N  =  26 

17  +  5IPTS 

1 

-  15, 

N  =  27 

OMEGA(N) 

See  card 

Floating 

16 

-  30, 

N  -  28 

12  +  5IPTS 

El  5.  7 

31 

-  45, 

N  =  29 

3* 

1 

-  5 

N  PRINT 

If  N PRINT  =  0 

Integer 

initial  data  plane 

15 

is  not  printed  out 

6 

-  10 

MODPLN 

A  plane  of  output 

Integer 

data  is  printed  out 

15 

every  MODPLN 
integration  steps. 

11-15  NPLOUT  A  plane  of  output  Integer 

data  is  printed  out  15 

on  binary  tape  4 
every  NPLUT 
integration  steps. 

16  -  20  NITER  Maximum  number  Integer 

of  cycles  allowed  15 

for  internal 
iteration  loops. 

(usually  NITER  -  10) 

21-25  NPLMAX  Maximum  integration  Integer 

plane  number  to  be  15 

calculated.  If  this  is 
a  restart  of  a  partially 
completed  case,  this 
number  should  be 
equal  to  number  of 
steps  already  completed 
+  number  of  steps 
desi red. 
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TABLE  A-I  (Cont.  ) 
TABLE  OF  INPUT  DATA  CARDS 


Card 

Number 

Column 

Number 

MNEMONIC 

Item  Entered  R 

26  -  30 

IPTMAX 

Maximum  number 
of  rings  allowed 

15  -  IPTMAX  -  21 

1  -  15 

CSST 

CSST  controls  the 
addition  of  a  new 
ring  of  streamlines 
at  the  shock.  A 
value  close  to  1.0 
is  customary  but 
this  may  be  in¬ 
creased  if  less 
detail  is  requ;  red 
near  the  shock. 

16  -  30 

CONST  1 

Step  size  control. 

The  step  size  is 
multiplied  by 

CONST  1  each  time 
it  is  computed. 

Usually 

0.  6  —  CONST  1  £1.0 

31  -  45 

CONST2 

Step  size  control. 
Limits  step  size  when 
a  new  ring  of  field 
points  are  added  at 
shock.  Usually 
CONST2  1.5 

46  -  60 

ERROR 

See  card  9.  columns 
16-30 

60  -  75 

XMAX 

Maximum 
x-coordinate  to  be 
computed. 

5=:* 

1  -  5 

ISW  1 

Control  constant  for 
internal  dumps.  See 
Section  A.  3  for  sub¬ 
routines  involved. 
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Number 

presentation 

Integer 

Ic 

F  loating 
E15.7 


F  loating 
E  1  5.  7 


F  Joating 
El  5.  7 


F  loating 
El  5.  7 

F  loating 
El  5.  7 

Integc r 
IS 


TABLE  A-I  (Cont.  ) 
TABLE  OF  INPUT  DATA  CARDS 


Ca  rd 
Number 

Column 

N  um  b  e  r 

MNEMONIC 

Item  Entered  Re] 

Number 

or  es  entation 

6  - 

10 

ISW  2 

Same  as  ISW  1 

Intege  r 

1  5 

1  1  - 

1  5 

ISW  3 

Same  as  ISW  1 

Intege  r 

15 

16  - 

20 

ISW  4 

Same  as  ISW  1 

Intege  r 

1  5 

21  - 

25 

ISW  5 

Same  as  ISW  1 

Intege  r 

1  5 

26  - 

30 

ISW  6 

Same  as  ISW  1 

Integer 

15 

6* 

1  - 

1  5 

RBDY 

Spherical  nose 
radius  in  cm  or  ft 

Floating 

E15.7 

16  - 

30 

SLOPE 1 

Slope  of  conical 
afterbody  for 
circular  cone 
sphere  combination 
(degrees) 

Floating 

E15.7 

31  - 

45 

SLOPE2 

Slope  of  expansion 
side  for  elliptical 
cone  at  0°  angle  of 
attack  (degrees) 

Floating 

El  5.  7 

46  - 

60 

SLOPES 

Slope  of  compression 
side  for  elliptical 
cone  at  0°  angle  of 
attac  k  (degrees) 

Floating 

El  5.  7 

61  - 

75 

BDYTYP 

Define  body  type  = 

1  .  0  for  a  spherical 

Floating 

E  1  5.  7 

nose  and  right 
circular  cone 
afterbody 

=  LO  for  a  spherical 
nose  and  an  elliptical 
afterbody  with  major 
axis  along  £  -direction 

=  3.0  for  a  spherical 
nose  and  an  elliptical 
afterbody  with  minor 
axis  along  z-direction 
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uiiiiiiiiimiiiiiiiiiiiiiiiifiiii 


Figure  A  3  INPUT  FORMAT  (Cont) 


This  completes  the  description  of  the  input  deck.  For  a  new  case,  a  total 
of  2  1  +  5IPTS  cards  are  required  for  input,  where  IPTS  is  the  number  of  data 
points  on  the  initial  data  ray.  For  a  restart  of  a  previous  problem,  a  total 
number  of  6  data  cards  are  required. 
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A.  3 


OUTPUT 


The  first  page  of  output  from  the  program  begins  with  the  run 
number  centered  under  the  heading.  Following  this  come  in  order:  one 
line  of  comment  (HEAD),  the  option  controls,  such  as  NSPEC,  MODEL.  , 
etc.  ,  the  free  stream  conditions,  the  free  stream  compositions  and  the 
input  body  description  data.  This  first  page  is  included  initially,  and  also 
if  the  present  run  is  a  restart  of  a  previously  incomplete  case. 

Following  the  initial  page  of  output  the  values  of  the  point  data 
for  a  plane  of  output  is  usually  printed  next.  The  plane  printed  is  the  initial 
plane  if  the  input  control  constant  N  PRINT  ■  1.  Otherwise,  every  MODPLN 
planes  are  printed  out,  where  MODPLN  is  an  input  constant.  The  output 
is  arranged  by  rings  starting  at  the  body  and  ending  with  the  points  on 
the  shock  surface  in  the  present  plane.  The  form  of  the  output  for  each 
point  in  the  plane  is  shown  in  Table  A-II.  The  printout  for  each  plane  is 
preceded  by  a  heading  (written  by  subroutine  OUTPUT)  labeling  the  form 
of  the  following  data.  Note  that  the  points  with  the  same  index  I  are  said 
to  lie  on  a  "ring"  and  points  with  the  same  index  N  are  said  to  lie  on  a 
"ray".  Although  these  were  circular  rings  and  straight  line  rays  in  the 
initial  plane,  they  distort  somewhat  as  the  computations  proceed  down¬ 
stream.  However  the  names  ring  and  ray  are  maintained  for  simplicity 
here. 

Other  output  can  be  obtained  if  any  or  all  of  the  input  control 

constants  ISW1 ,  ISW2, _ ,  ISW6  are  equal  to  1.  In  this  case  internal 

output  statements  print  out  certain  data  within  the  various  subroutines. 

These  are  helpful  in  determining  the  cause  of  any  problems  which  may 
develop  during  a  run.  While  the  forms  of  these  optional  outputs  will  not 
be  given  here,  they  are  easily  obtained  by  reference  to  the  aapropriate 
FORTRAN  listing.  Each  optional  output  is  preceded  by  the  subroutine 
name  of  the  subroutine  from  which  it  was  written.  Following  is  a  list  of 
the  subroutines  which  have  internal  dumps  and  the  control  constants  which 
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initiate  their  use.  Caution  should  be  exercised  when  using  the  internal 
dumps  as  they  generate  large  amounts  of  output  in  a  comparatively  short 
time. 


1. 

ISW1  =  1 

Controls 

2. 

ISW2  =  1 

Controls 

3. 

1SW3  =  1 

Controls 

4. 

ISW4  =  1 

Controls 

E>. 

IS  W5  =  1 

Controls 

6. 

ISW6  =  1 

Controls 

MCMAIN,  ROTATE 
POINT 

FIELD,  SECT,  POINT  B 
CDELTA.  FLOW,  SETUP,  SHOCK 
(not  used) 

BODY,  SECT  B 


Table  A-II 

FORM  OF  OUTPUT  FOR  EACH  PLANE 


RESULTS  FOR  N  PLANE  =  (N  PLANE) 


COLUMN  NUMBER 
4  5 


MACH 

NUMBER 


(DEGREES) 


ALPHA  1  ALPHA  2  ALPHA  3 


’NOTE:  LINE  4  APPEARS  0NLv  WHEN  THE  POINT  BEING  OUTPUT  IS  ON  THE  SHOCK 
(i.e.  THE  MAXIMUM  RING  NUMBER) 
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A.  4 


DATA  STORE 


The  data  for  a  plane  of  points  are  stored  in  a  large  data  matrix 
A1  (25,  30,  1  5)  which  allows  for  up  to  25  rings  of  points  with  30  points 
(rays)  per  ring  and  15  functions  per  point.  For  program  simplicity,  successive 
blocks  of  7500  storage  locations  have  been  equivalenced  to  a  named  function 
which  in  turn  is  used  throughout  the  program  in  place  of  the  large  matrix  A1 
whenever  a  certain  one  of  the  functions  stored  in  A1  is  required.  Following 
is  a  list  of  these  new  MNEMONICS  giving  their  position  in  A1  and  their 
significance 

Starting 

Function _ Nmemonic _ (Position  in  Al) _ Definition _ 


1 

VY(25,  30) 

Al(l, 1,  1) 

Y -coordinate 

2 

VZ(25,  30) 

Al  (1 ,  1 , 2) 

Z -coordinate 

3 

VP(25,  30) 

Al  (1 ,  1 ,  3) 

pressure 

4 

VTHETA(25,  30)  = 

Al(l,  1,4) 

flow  angle,  See  Fig.  1 

5 

VPSI(25,  30) 

Al(l,  1,  5) 

flow  angle,  See  Fig.  1 

6 

VQ(25,  30) 

Al(l,  1,6) 

velocity 

7 

VBETA(25,  30) 

Al(l,  1,  7) 

Mach  angle 

8 

VF(25,  30) 

Al(l,  1,8) 

not  used  for  frozen  solution 

9 

VGAMMA(25,  30)  = 

Al(l,  1,9) 

specific  heat  ratio 

10 

VMW(25,  30) 

Al(l,  1,  10) 

Molecular  Weight 

1  1 

VRHO(  25,  30) 

Al(l,  1,11) 

Density 

.  2 

VDSTL(25,  30) 

Al  ( 1 ,  1,12) 

distance  along  streamline 

13 

VX(25,  30) 

Al(l,  1,  13) 

X-coordinate 

14 

VR( 25,  30) 

A 1  ( 1 ,  1,14) 

Polar  RADIUS 

15 

VOMEGA(25,  30)  = 

Al(l,  1,  15) 

Polar  Angle 

A.  4.  1 

Program  Stops 

When  an  error  (other  than  Fortran)  occurs  within  the  program 
the  program  stops  with  an  appropriate  error  message  and  full  octal  core 
dump.  The  error  message  refers  back  to  the  statement  number  of  the 
statement  and  the  subroutine  in  which  the  error  occurs.  In  most  cases  the 
error  is  self  explanatory. 
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A.  5 


DISCUSSION  OF  INDIVIDUAL  SUBROUTINES 


A  detailed  synopsis  of  the  main  program  and  the  individual 
subroutines  is  presented  in  the  present  section  in  order  to  point  out 
techniques  and  procedures  which  have  not  previously  been  discussed 
elsewhere  in  the  report. 

A.  5.  1  M  CM  AIN 

The  main  program  controls  the  overall  logic  of  the  whole 
integration  scheme,  by  first  locating  planes  of  symmetry  and  then  shifting 
from  one  point  to  another  throughout  the  initial  data  plane  in  order  to  compute 
new  data  points  in  the  next  plane.  The  data  is  arranged  so  that  N  data 
points  lie  along  I  rings  from  the  body  to  the  shock,  N  being  constant  for 
all  I.  Starting  on  the  expansion  side  of  the  body  at  a  point  located  on  the 
shock  (the  Ith  ring)  in  the  plane  of  symmetry  (  A/  =  l)A/new  shock  points  in 
the  next  data  plane  are  successively  computed.  Control  is  then  shifted  to 
the  ring  of  body  points.  After  N  new  body  points  are  computed,  rings  of 
field  points  between  the  body  and  the  shock  are  computed  sequentially.  As 
the  calculations  progress,  data  points  on  the  new  data  plane  are  written 
over  those  in  the  initial  plane  which  are  no  longer  needed.  In  this  way,  as 
much  storage  space  as  possible  is  saved,  since  no  more  than  one  surface 
need  be  stored  at  a  time. 
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When  the  ring  of  shock  points  is  again  reached,  two  options  are 
possible.  If  a  control  parameter  called  IADD  is  not  O,  a  new  ring  of  field 
points  corresponding  to  the  streamlines  through  the  shock  points  is  added 
in  the  new  data  plane,  otherwise  the  computation  of  new  points  is  stopped. 

In  this  way,  as  the  shock  layer  grows  along  the  body,  the  mesh  size 
through  the  layer  can  be  maintained  fairly  uniform  without  resorting  to 
interpolation  procedures  in  order  to  add  new  points  to  the  field.  Figure  27 
shows  some  typical  streamlines  that  would  be  added  at  the  shock  in  a  plane 
of  symmetry.  If  a  new  ring  of  field  points  is  added  at  the  shock,  the  ring 
of  shock  points  in  the  new  data  plane  is  indexed  up  to  I  +  1  and  this  data 
plane  then  becomes  a  new  initial  plane  from  which  another  plane  can  be 
calculated.  New  planes  are  calculated  in  this  manner,  until  either  a 
maximum  number  of  planes  have  been  determined,  or  a  maximum  * 
coordinate  is  reached. 

The  main  program  also  has  the  option  of  dropping  whole  rings 
of  data  points  whenever  desired.  Presently  this  is  accomplished  by  first 
testing  whether  the  total  number  of  rings  J  is  equal  to  ihe  maximum 
number  allowable  IPTMAX,  an  input  constant  to  the  program.  If  J  equals 
IPTMAX,  and  the  mesh  spacing  near  the  shock  is  such  that  a  new  ring  of 
shock  points  would  be  added  in  tha  next  data  plane  (agai.  in  order  to  keep 
mesh  spacings  fairly  "niform)  then,  every  other  ring  of  field  points  between 
the  body  and  the  shock  is  dropped  from  the  computation.  This  technique 
takes  advantage  of  the  fact  that  since  the  calculations  proceed  along  the 
x-axis  and  the  flow  gradients  decrease,  less  detail  is  required  to  completely 
specify  the  flow  field.  Hence,  the  overall  computation  procedure  can  be 
drastically  speeded  up  --  both  from  the  points  of  view  of  having  now  fewer 
points  to  calculate,  and  having  also  larger  mesh  widths  in  the  initial  data 
plane,  which  in  turn  allow  for  larger  axial  step  sizes  AX  .  For  many 
problems  it  may  be  desirable  to  drop  rings  other  than  these  throughout 
the  flow  field.  This  is  especially  true  for  some  bodies  at  high  angles-of- 
attack  when  the  spacings  between  rings  near  the  bocy  on  the  compression 
side  may  tend  to  be  small.  By  checking  spacings  between  the  rings  in  say 
the  plane  of  symmetry  on  the  compression  side  of  the  body,  any  specified 
rings  can  be  dropped  from  the  procedure. 
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MCMAIN  calls  subroutines  OUTPUT,  STOP,  STEP,  FIT, 
BODY,  FIELD  and  SHOCK  in  that  order.  The  function  of  each  of  these 
routines  is  given  below.  Since  none  of  the  actual  computations  are  per¬ 
formed  in  MCMAIN,  only  slight  changes  to  account  for  any  extra  functions 
which  are  added  will  be  required  in  order  to  include  different  thermo¬ 
chemical  models  in  the  program. 

A.  5.2  INPUT 

This  subroutine  reads  all  the  input  data  that  is  required  by  the 
program  from  a  tape  where  it  had  previously  been  written  in  the  proper 
order  and  format.  INPUT  has  essentially  two  separate  functions, 
depending  upon  what  data  is  read  from  the  tape.  If  the  data  is  that 
obtained  from  a  blunt  body  solution  ale  lg  one  line  only,  it  is  first  non- 
dimensionalized  using  the  methods  given  in  Paragraph  A.  3,  and  then  trans¬ 
ferred  to  several  different  azimuthal  planes  using  the  fact  that  the  initial 
data  should  be  axisymmetric.  In  this  case,  the  input  data  is  also  written 
on  BCD  output  tape  for  checking  purposes. 

If  the  data  is  that  read  from  a  binary  tape  (auxiliary  tape  2) 
which  has  previously  been  w ritten  by  the  program  (see  subroutine  STOP) 
it  then  represents  a  continuation  of  a  previously  run  problem  and  is 
already  in  the  form  necessary  for  use  by  the  program.  Hence,  control 
is  transferred  back  to  MCMAIN  immediately. 

The  initial  point  data  thus  acquired  by  INPUT  is  stored  in  a 
common  three-dimensional  array  labeled  Al.  Any  change  in  the  thermo¬ 
chemical  model  would  be  reflected  in  a  corresponding  increase  in  the  size 
of  Al  in  storage.  The  correspondence  between  the  elements  of  Al  and 
the  nomenclature  used  in  the  program  is  given  in  Paragraph  A.  4.  Calls 
SMALLH  and  PAUSE. 

A.  5.  3  STEP 

The  STEP  size  A Tt>  between  successive  data  planes  normal  to 
the  x-axis  is  a  function  of  the  spacing  of  the  points  on  the  base  plane,  the 
magnitude  and  direction  of  the  velocity  field,  and  the  speed  of  sound 
throughout  the  field.  Here,  A*  is  chosen  so  that  the  Mach  forecone  from 
any  new  point  intersects  the  base  plane  in  such  a  manner  that  its  domain 
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of  dependence  does  not  contain  any  initial  points  adjacent  to  the  original 
bare  point.  This  can  in  general  be  assured  by  taking 

Ax  <  D/ST  (cosip,^* 

where  DIST  is  the  minimum  mesh  point  spacing,  Am3iX  is  t*ie  maximum 
Mach  angle  and  cf>maK  is  t*ie  maximum  angle  which  the  velocity  vector 
makes  with  the  x-axis.  Usually,  however,  it  will  not  be  necessary  to 
check  every  point  in  the  initial  plane  in  order  to  satisfy  (A-4).  In  the  present, 
version  of  STEP,  only  points  near  the  body  and  the  shock  in  the  geometrical 
planes  of  symmetry  are  checked  in  order  to  determine  by,  .  At  as 
obtained  from  (A-4)  is  then  multiplied  by  an  input  constant  (CONST)  to 
yield  the  final  va.ue  of  At  .  For  problems  where  the  maximum  values 
of  the  Mach  angle  and  flow  angle  are  expected  off  the  planes  of  symmetry 
(e.g.  ,  an  ellipsoidal  body  at  angle  of  attack),  CONST1  can  be  appropriately 
adjusted  so  that  the  computation  can  proceed. 

Another  function  of  STEP  is  to  test  the  spacing  of  data  points 
near  the  shock  surface  to  determine  whether  or  not  a  new  str<  amline  will 
be  added  at  the  shock.  If  new  streamlines  are  required,  the  test  constant 
IADD  (see  MAIN)  is  set  equal  to  one  --  otherwise,  IADD  equals  zero. 

A.  5.  4  FIT  And  LSTSQ 

Interpolating  surface  fits  through  nine  initial  data  points  are 
obtained  in  this  subroutine  along  the  lines  discussed  in  Step  1  of  Paragraph 
3.3.2.  The  functions  fit  are  the  pressure,  velocity  and  two  flow  angles. 
Since  the  flow  angles  9  and  f  vary  considerably  around  a  given  body 
(e.g.,  at  any  x-locaticn  for  an  axisymmetnc  flow,  9  varies  from  a 
maximum  along  the  i^-axis,  to  zero  along  the  ^-axis)  a  transformation 
to  new  spherical  angles  uj ,  x  is  first  performed.  The  justification  for 
for  this  lies  somewhat  in  the  fact  that  for  many  three-dimensional  problems 
considered,  the  flow  angle  u>  will  vary  mainly  with  one  of  the  independent 
variables,  while  x>  will  vary  mainly  with  the  other.  The  relationship 
between  the  angles  a/,  x  and  the  angles  6  ,  is  shown  in  Figure  ( A-4). 
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The  independent  variables  presently  used  in  the  fitting  technique 
are  the  polar  coordinates  (  r  ,J1)  in  the  yj.-plane.  Here, 


n. 


(A-5) 


is  the  angle  between  the  y-axis  and  the  point.  While  these  coordinates 
are  generally  best  for  most  "rounded"  bodies, (e.g.  ,  axisymmetric  or 
ellipses  of  moderate  eccentricities),for  bodies  with  locally  two  dimensional 
regions  (e.g.  ,  a  "flat"  ellipse)  perhaps  a  combination  of  cartesian 
coordinates  (y,  ^ )  over  the  flat  portion  and  polar  coordinates  (ntJ~L)  over 
the  rounded  portions  would  be  better.  In  this  case,  the  focal  point  of  the 
ellipse  would  be  chosen  as  the  center  of  the  (y,  j.)  coordinate  system. 

This  can  better  be  determined  after  gaining  more  experience  with  the 
program.  Call  HETA,  SIMSOL. 

A.  5.  5  FIELD  (IP,NP) 

Given  an  initial  field  point  (IP,  NP),  the  step  size  Ax  and  the 
results  from  subroutine  FIT,  the  FIELD  subroutine  calculates  the 
coordinates  and  flow  properties  at  a  r.ew  field  point.  To  do  this,  subroutine 
POINT,  STREAM,  SECT,  INTPLT  and  SETUP  are  used  throughout  an 
iteration  loop  to  converge  the  solution  to  second  order  in  the  step  size, 
along  the  lines  proposed  in  Paragraph  3.  3.  2. 

As  an  error  check  in  the  modified  Euler  iteration  scheme,  a 
limit  of  NITER  cycles  is  set  on  the  looping.  If  the  convergence  test  on  pressure 


p<» > 

p(n) 


<  0.00001 


(A-6) 


is  not  satisfied,  within  these  NITER  cycles,  an  error  statement  is  written 
off-line  and  control  is  returned  to  the  main  program.  Calls  SIMSOL,  INTPLT, 
POINT,  STREAM,  POINTB,  SCKFN,  SHIFT,  SECT,  SETUP. 

A.  5.6  BODY  (IP,  NP) 

Given  an  initial  body  point  (IP,  NP),  the  BODY  subroutine 
calculates  the  flow  properties  at  a  point  where  the  body  streamline  through 
(1P,NP)  meets  the  next  data  plane.  To  accomplish  this,  BODY  calls 
exactly  the  same  subroutines  as  FIELD  plus  two  more  given  the  mnemonics 
BODYFN  and  SECTB.  The  procedure  is  essentially  that  given  in  Paragraph 
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3.3.3,  where  now,  two  iteration  loops  are  required  in  order  to  converge 
the  solution  to  second  order  in  the  step  size. 

The  inner  loop  uses  Newton  iteration  to  converge  the  flow  angles 
6  and  (/>  such  that  the  boundary  condition  of  flow  tangency  at  the  surface  is 
solved  with  the  absolute  limits 

0<n)- 0(n'n  <  fO'a  (A- 7) 

and 

yto-f)  <  l0-8  (A-8) 

If  either  of  the  Equations  (A-7)  and  (A-8)  is  not  satisfied  within  20  cycles  an 
error  statement  is  again  printed  off-line  and  the  computations  allowed  to 
continue. 

The  outer  iteration  is  again  a  modified  Euler  type  on  pressure 
and  contains  the  same  error  checks  and  convergence  criterion  as  does 
subroutine  FIELD. 

A.  5.  7  SHOCK  (IP,  NP) 

Given  an  initial  shock  point  (IP,NP),  this  subroutine  calculates 
the  coordinates,  flow  properties  and  direction  cosines  of  the  shock  normal 
at  a  new  point  on  the  shock  surface  in  the  next  data  plane.  The  method  is 
that  given  in  Paragraph  3.  3.4.  In  the  process,  SHOCK  uses  subroutines  FLOW, 
CDELTA  and  SETUP,  whose  various  functions  are  explained  below,  to 
converge  the  solution  to  again,  second  order  in  the  step  size.  Here  again, 
the  convergence  tests  on  pressure  are  the  same  as  those  used  in  FIELD 
and  BODY  with  the  appropriate  error  statement  included.  Calls  SHIFT, 

PAUSE,  POINT. 

A.  5.  8  INTPLT 

Given  the  cartesian  coordinates  (y ,  j )  of  a  point  P  in  the  base 
plane,  INTPLT  calculates  its  polar  coordinates  (t',,n)  (See  Equation  (A--1)) 
and  uses  these  in  the  equations  defined  by  FIT  in  order  to  determine 
pressure,  velocity  and  llow  angles  B  and  f  at  P.  &  and  are  obtained 
from  uJ  and  £  (the  actual  angles  fit)  by  using  the  inverse  transformations 
(See  Figure  (A-4). 
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and 


=  ta.n~f(6anx 5  sin  u)) 


(A-9) 


£  s  ifa-o  *  '  f ^  )  (A- 10) 

CALLS  HETA 
A.  5.9  STREAM 

The  compatibility  equations  written  along  the  streamline  are 
solved  in  STREAM  using  in  this  case,  modified  Euler  integration.  For 
equilibrium  and  frozen  flow,  this  same  integration  technique  can  be  used. 
However,  for  nonequilibrium  flow  a  more  accurate  technique  using  the 
methods  discussed  in  Paragraph  3.  3.  2  would  require  an  entirely  new 
STREAM  routine.  Calls  PAUSE. 
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A.  5.  10  SECT 


Given  the  flow  properties  at  the  new  point  and  the  bicharacteristic 
parametric  angle  S  ,  this  subroutine  solves  for  the  coordinates  of  the  new 
base  point  in  the  initial  plane  in  the  manner  indicated  in  Step  3  of  Paragraph 
3.3.?.. 

A.  5.  11  SETUP 

In  this  subroutine  the  coefficients  of  the  compatibility  equations 
written  along  the  bicharacteristics  are  obtained  in  the  manner  indicated 
in  Step  5  of  Paragraph  3.  3.  2.  For  nonequilibrium  flow,  the  only  change 
required  here  would  be  to  add  the  term  -sm  A Lt  to  the  calcul¬ 

ations  along  the  bicharacteristic  i  .  Calls  META. 

A.  5.  1  2  SECTB 

In  SECTB,  the  new  body  point  is  located  as  the  intersection  of 
the  plane  containing  the  normal  to  the  body  and  the  streamline  through  the 
initial  body  point  with  the  equation  of  the  body  surface  and  the  next  plane. 
Since  the  equation  of  the  body  surface  usually  changes  for  every  problem, 
the  equations  are  solved  using  a  second  order  Newton  process.  Calls 
BDYFN,  PAUSE. 

A.  5.  13  BODYFN 

The  equation  of  the  body 

5fx,  y,  p)  -  O 

and  the  direction  numbers  of  the  body  normals  are  supplied  by  BODYFN, 
which  must  be  rewritten  if  a  sphere  cone  is  not  desired.  Here  the  conical 
afterbody  maybe  circular  or  elliptical  as  specified  on  input.  The  body 
functions  are  written  with  respect  to  a  body  centered  coordinate  system 
with  the  X-axis  as  the  body  axis.  Eventually,  BODYFN  should  be  replaced 
by  a  general  routine  which  allows  for  either  an  analytical  description  or  a 
numerical  surface  fitting  technique  to  describe  the  body.. 

A.  5.  14  SCKFN 

Three  adjacent  shock  points  on  the  shock  surface  in  the  initial 
plane  are  fit  with  an  exact  quadratic  fit  in  the  polar  coordinates  (  r,n)  by 
SCKFN.  Calls  HETA. 
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A.  5.  15  CDELTA 


This  subroutine  locates  three  bicharacteristics  passing  back 
from  the  new  shock  point  to  the  initial  plane,  using  the  method  presented 
in  Step  4  of  Paragraph  3.  3,  4,  A  Newton-Raphsen  iteration  scheme  is  used  to 

solve  for  the  base  coordinates  of  each  of  two  bicharacteristics  imbeded  in 
the  shock  surface.  Then,  a  third  bicharacteristic  whose  direction  is 
represented  by  the  average  of  these  wo  is  located,  CDELTA  also  calls 
INTPLT  in  order  to  determine  the  base  point  flow  properties.  Calls 
SCKFN,  POINT B,  SIMSOL,  HETA,  INTPLT. 

A.  5.16  STOP 

Subroutine  STOP  has  two  basic  functions.  First,  a  logical 
binary  record  containing  all  the  data  which  would  be  necessary  to  restart 
the  computations  is  written  for  every  //planes,  where  hj  is  an  input 
constant  to  the  program.  Then,  if  either  a  given  maximum  number  of 
planes  have  been  computed,  or  a  given  maximum  X.- coordinate  has  been 
reached,  the  computational  procedure  is  halted,  otherwise,  control  is 
transferred  back  to  MCMAIN. 

A.  5.17  OUTPUT 

Output  data  from  the  program  is  written  for  every  A/  planes, 
where  here,  /V  is  the  same  input  constant  which  was  used  in  STOP. 

OUTPUT  can  easily  be  changed  to  provide  any  output  quantity  desired, 
without  affecting  the  reaminder  of  the  program.  The  form  of  the  output 
obtained  at  each  point  in  the  plane  in  the  present  version  of  OUTPUT,  is 
described  in  Paragraph  A.  3.  OUTPUT  calls  POINT  which  prints  out  the 
point  data  for  each  point  in  the  plane. 
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A.  5.  1  8  SIMSOL 

Any  set  of  yi  linear  inhomogeneous  equations  in  unknowns 
can  be  solved  in  SIMSOL.  If  the  determinant  of  the  coefficients  vanishes, 
an  appropriate  error  statement  is  printed  out. 

A.  5.19  DIVIDE 

DIVIDE  is  used  to  drop  every  other  ring  of  data  points  from  the 
shock  to  the  body. 

A.  5.  20  ROTATE 

ROTATE  is  used  to  rotate  the  coordinates  such  that  all  the  data 
points  lie  in  a  plane  x  =  constant.  The  option  which  utilizes  this  subroutine 
is  not  presently  working  properly.  Calls  HETA, 

A.  5.21  POINT 

POINT  is  used  to  write  out  the  data  at  a  given  data  point.  The 
form  of  this  output  is  given  in  Paragraph  A.  3. 

A.  5.  22  POINTB 

POINTB  is  used  to  provide  extra  output  data  at  the  base  of  each 
bicharacteristic.  It  is  utilized  mainly  for  debugging  purposes. 

A.  5.23  PAUSE 

If  an  error  occurs  which  would  make  it  impossible  to  continue 
the  procedure,  an  appropriate  error  statement  is  written  by  PAUSE,  and 
the  program  is  terminated. 

A.  5.  24  HETA 

Determine  an  angle  6  ,  in  radians  {-  -  0  ^  -j  )  when  given 

the  opposite  and  adjacent  sides  of  a  right  triangle. 
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A.  5.  25 


FLOW 


FLOW  is  used  to  compute  the  conditions  behind  an  oblique  shock 
given  t lie  shock  angle,  ,  and  the  tree-stream  conditions  at  each  point. 

For  all  thermodynamic  options  the  computations  utilize  a  Newton-Raphsor. 
procedure  to  match  the  total  enthalpy  behind  the  shock  with  that  of  the  free 
stream.  FLOW  calls  SMALLH  and  PAUSE, 

A.  5.  26  SMALLH 

SMALLH  is  presently  used  to  provide  harmonic  oscillator  approx¬ 
imations  to  the  species  thermodynamic  functions  enthalpy,  A>  ,  specific  heat, 

y/ 

Cf>j,  and  free  energy^.  .  The  harmonic  oscillator  data  is  contained  within 
the  program  in  the  form  of  data  statements  with  a  maximum  of  8  species 
allowed.  These  are  presently,  O,  N,  E  ,  A,  O^,  N^,  NO,  NO  +  ,  in  that 
order.  One  data  card  is  read  in  INPUT  in  order  to  choose  the  proper  data 
for  the  thermochemical  system  being  run.  A  description  of  this  card  is 
included  in  the  section  on  the  discussion  of  the  program  input  data. 

A.  5.  27  STMN 

Control  program  which  is  used  to  read  run  label  and  control  the 
overall  execution  of  the  program.  STMN  calls  INPUT  and  MCMAIN. 

A.  5.  28  SHIFT 

SHIFT  is  used  to  shift  the  point  data  from  location  II  to 
location  12. 

A.  5.  29  COMMON  DECK 

Must  be  used  with  all  subroutines  except  SIMSOL,  HETA  and 

PAUSE. 
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SYMBOL  DEFINITION 


A.  6 

Following  is  a  list  of  symbols  which  are  used  throughout  the  program. 
Arrangement  is  in  the  order  the  variable  appears  for  each  labeled  common 
block  followed  by  the  list  of  special  variables  used  within  each  subroutine. 

The  asterisk  designates  input  quantities. 


Mnemonic  Variation 

Symbol 

Comment 

COMMN/  FREE 

:':P0  ^  PRES 

/ 

'P'co 

2 

Free-stream  pressure  in  dynes/cm 

or  lb/ ft^ 

RHOO  -  DENS 

Pi 

3 

Free-stream  density  in  grams/cm 

or  slugs  / ft.3 

"IK)  -  VEL 

ui 

Free-stream  velocity  in  cm/sec  or 

ft  /  sec 

1  WMW  0  -  WGHT 

MWj 

Free-stream  molecular  weight  in 

grams  lbs. 

mole  mole 

'  TO  -  TEMP 

r' 

Free-stream  temperature  in  “K 

"HINF 

Free-stream  total  enthalpy  normalized 

by  KTJ/mw; 

TINF 

T. 

Free-stream  temperature  normalized 

by  Ti 

PAM 

P- 

Free-stream  pressure  normalized 

by  pU  u‘m% 

RHAM 

P» 

Free-stream  density  normalized  by 

K 

UAM 

Free-stream  velocity  normalized  by 

U  ' 

OO 

WTAM 

Free-stream  molecular  weight 
normalized  by  MW# 

"DIML 

Reference  length  in  cm  or  ft 

:  R  0 

—  / 

R 

Universal  Gas  constant  - -.^S  .or  — '  l*3  • 

moleA  mole- * 

BIGAM 

t 

90 

Free-stream  specific  heat  ratio 
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Mnemonic  Variation 

XLMBDA  TSCALE 

"'EPSAM(K) 

-  EPSINF(K) 

GAMAM(J) 

-GAMINF(J) 

COMMON/THERMO 
HJO  (J) 

H(J) 

F  J  (J) 

CP(J) 

EPSEQ(K) 

NATOMS(J) 

COMMON  /  DATA3D 
A  1(1 ,  J ,  K) 

ALF2(N) 

ALF3(N) 

ALFL2(N) 

ALFL3(N) 

‘GAMMA I  (I,  J) 

EPS  I  (I,  K) 


Symbol  Comment 

t  » 

,  MW.  Ul, 

A  Defined  as  .  /  * — “ — 

*  V 

Free-stream  vibrational  energy  for 
kth  vibrator  normalized  by 

f  Free-stream  concentration  of  Jth 

o»j 

species  normalized  by  2_  t 

<30  ‘ 

J  J 


Energy  of  formation  of  jth  species 
nondimensionalized  by  R  'TJ  / wwj, 


h  ■ 

j 


Enthalpy  of  jth  species  nondimension¬ 
alized  by  R'tJ 


Free-energy  of  jth  species  non¬ 
dimensionalized  by  R't^,/wwJ 

Specific  heat  of  jth  species  non¬ 
dimensionalized  by  r'/m wj 


Equilibrium  vibrational  energy  of 
kth  vibrator  normalized  by  R'T^/bWJ^ 


NJ 


Number  of  atoms  in  species  j. 


Storage  matrix  for  point  data 
I  =  ring  number,!  ray  number, 

K  -  variable  number 

Direction  cosine  of  shock  normal  for 
shock  point  (N)  along  y-axis 

Direction  cosine  of  shock  normal  for 
shock  point  (N)  along  z-axis 

ALF2(N)  for  previous  shock  point  (N) 

ALF3(N)  for  previous  shock  point  (N) 


Concentration  for  ring  I  and  species  J 


Vibrational  energy  for  ring  I  and 
vibrator  K 
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Mnemonic  Variation  Symbol 
COMMON /BOSS  3D 


Comment 


AM{9,  10) 

ALPHAL 

ALPHAN 

BETAMX 

CCP 

CSIGMA 

SSIGMA 

CPSIFS 

SPSIFS 

CTHEFS 

STHEFS 

'cONSTl 


Matrix  used  for  storing  coefficients 
of  linear  equations  to  be  solved 
simultaneously 

Angle  between  present  x-axis  and 
body  axis  when  plane  rotation  is  used 

Angle  between  normal  to  rotated 
plane  and  body  axis 

Maximum  Mach  angle 

Specific  heat  ratio  (nondimensional) 

Cosine  of  shock  angle 

Sine  of  shock  angle 

Cosine  of  free  stream  flow  angle 

Sine  of  free  stream  flow  angle 

Cosine  of  free  stream  flow  angle 

Sine  of  free  stream  flow  angle 

Input  constant  used  to  regulate 
step  size, 


CONST2 


DALPHA 

DLMAX 

DXMIN 


Input  constant  used  to  decrease  step 
step  size  when  new  field  points  are 
added  at  the  shock 

Angle  of  rotation  of  data  plane  when 
coordinate  rotation  is  employed 

Not  used 

Minimum  step  size  allowed 


IPTS 


Number  of  data  rings  for  present  step 


ICONST(IO) 


Integer  matrix  (not  used) 


:  IPTMAX 


Maximum  number  of  data  rings 
allowed 
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Mnemonic  Variation  Symbol 

"iswi 
"lSW2 
*  ISW  3 
ISW4 
: :  IS  W  5 
ISW  6 

:  IT  A  PE 
*MDELAY 

:':MODPLN 
"MODEL 

"nc 

NEWRUN 

"NP 180 
"NPCMAX 

"nplout 
"nprint 

1 NROTAT 

"NPTS 
NPTMAX 
NPLANE 

"NR 


Comment 

Integer  test  constants  used  to  control 
internal  prograni  print  out 

Integer  Control  constant  (see  input) 

Integer  constant  always  set  =  3 
on  input 

See  description  in  input  section 
See  description  in  input  section 

See  description  in  input  section 

Control  constant  used  to  test  whether 
or  not  a  subroutine  has  been 
previou; ly  entered 

See  description  in  input  section 

See  description  in  input  section 

See  description  in  input  section 

See  description  in  input  section 

See  description  in  input  section 

See  description  in  input  section 
Not  used 

Number  of  current  integration  plane 
See  description  in  input  section 


1S4 


Mnemonic  Variation  Symbol 


Comment 


NRCORD 

'"NVIB 

NS  PEC 

PGRAD 
PP 
QQ 

'"SIGMA 

"SIGMA  1 
THETMN 
VARIBL(IO) 
xo 

YO 

'"XMAX 

"nospec(n) 

COMMON/AFT3D 
A2(J,  M,  N) 

AMUB(JJ) 

AMUT3 
AMUAV 

ALFT1  | 

ALFT2  f 
ALFT3  J 
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Not  used 

See  description  in  input  section 

See  description  in  input  section 

See  description  in  input  section 
Not  used 
Not  used 

Shock  angle 

Shock  angle  for  upstream  shock  point 

Minimum  flow  angle 

Real  marix  (not  used) 

x -coordinate  of  present  data  plane 

y -coordinate  of  last  shock  point  on 
present  data  plane 

Maximum  x-coordina^e  calculated 

See  description  in  input  section 

Matrix  used  in  setting  up  quadratic 
surface  fits. 

Mach  angle  at  bicharacteristic  (JJ) 
base  point. 

Mach  angle  at  new  point 

Average  Mach  angle  along  hicharacter- 
i  stic 

Temporary  storage  locations  for 
direction  cosines  to  shock  normal 


Mnemonic  Variation 

Symbol 

Comment 

AT32(K) 

Flow  angle,  0  ,  as  computed  by 

the  (  K  )  set  of  bicharacteristics 

AT32C 

Average  pressure  at  new  shock 
point  or  average  flow  angle,  0  , 
at  new  field  point 

AT33(K) 

Flow  angle,  f  ,  as  computed  by  the 
(K)  set  of  bicharacteristics. 

AT33C 

Average  flow  angle,  f  ,  at  new 
point 

AT31(K) 

Pressure  as  computed  by  the  (K) 
set  of  bicharacteristics 

AT31C 

Average  pressure  at  new  point 

AT3T(I,K) 

Storage  matrix  for  pressure  (It=l) 
for  flow  angle  0  (1  =  2)  and  for  flow' 

angle  y*  (1=3)  as  computed  from  the 
(K)  set  of  bicharacteristic s 

AVE(J) 

AVE(J)  is  equal  to  value  of  function 
(J)  for  middle  point  of  nine  points 
fit. 

AMUBT(J.K) 

Mach  angle  for  bicharacteristic  (J) 
in  the  (K)  set  of  bicharacteristics 

ALF 1 

*•» 

Direction  cosine  of  normal  to  shock 
with  respect  to  x-axis 

BO 

Body  function  B(x,  y,  z)  =  0  solved 
at  point  (x,  y,  z)  gives  remainder  BO. 

CDEL 

Cosine  of  parametric  angle  8  at 
new  point 

SDE.L 

) 

Sine  of  parametric  angle  6  at  new 
point 

CDELB 

Cosine  of  parametric  angle  S  at 
bicharacteristic  base  point 

SDELB 

Sine  of  parametric  angle  S  at  bi¬ 
characteristic  base  point 

1  S  6 


Mnemonic  Variation 

Symbol 

C  om  me  nt 

CAMU 

Cosine  of  Mach  angle 

SAMU 

AOTO  C/l) 

Sine  of  Mach  angle 

CAMUB 

Cosine  of  Mach  angle  at  bicharacter¬ 
istic  base  point 

SAMUB 

</*g) 

Sine  of  Mach  angle  at  bicharacteristic 
base  point 

CTHET 

XXHU  ( 6  ) 

Cosine  of  flow  angle  9 

STHET 

A Zru  (6) 

Sine  of  flow  angle  & 

CTHETB 

yCAd/  6) 

Cosine  of  flow  angle  6  at  bicharacter¬ 
istic  base  point 

STHETB 

Aisns  (6^) 

Sine  of  flow  angle  Q  at  bicharacter¬ 
istic  base  point 

CTHETO 

Cosine  of  flow  angle  9a 

STHETO 

Aorv  f  6,  ) 

Sine  of  flow  angle  9, 

CPSIO 

AW  (%) 

Cosine  of  flow  angle 

SPSIO 

AUu(K) 

Sine  of  flow  angle  HJ„ 

C01 

Direction  number  of  body  slope  with 
respect  to  x-axis  at  initial  body  point 

coz 

Direction  number  of  body  slope  with 
respect  to  y-axis  at  initial  body  point 

C03 

Direction  number  of  body  slope  with 
respect  to  y-axis  at  initial  body  point 

CT1 

CT2  ^ 

CT3  J 

Direction  number  of  body  slope  with 
respect  to  x?  y-and  z-axis  respectively 
for  new  body  point 

COEFN(N,  J) 

Coefficients  used  in  interpolating 
surface  fits  for  function  J. 

DELTA(J) 

* 

Parametric  angle  5  for  bicharactc r- 
istic  J. 
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Mnemonic  Variation  S ymbol 


Comments 


DELTAB(J)  Sb 


DELTBT(J,  K) 


DERFN(M,  N) 

DADX  ^ 
DADY  J- 
DADZ  J 

DADN 


DBDXAV 
DBDYAV  \ 
DBDZAV  J 

DALF2 


DALF3 


DELTAV 


DL(  J) 
DXNEXT 


DXDN  ■) 
DYDN  V 

DZDN  J 

El 
E2 
E3 
FI 
F  2 
F3 
G1 
G2 
G3 


Parametric  angle  (<5'b  )  at  base 
point  of  bicharacteristic  J. 

Parametric  angle  for  bicharacter¬ 
istic  (J)  in  set  of  bicharacteristics 
(K). 

not  used 

Derivative  of  speed  of  sound  with 
respect  to  x,  y  and  z  axes  respec¬ 
tively 

Derivative  of  speed  of  sound  witii 
respect  to  N-direction 

Average  derivative  of  body  function 
with  respect  to  the  x,  y  and  z  axes 
re  spectively 

Change  in  direction  cosine  of  shock 
normal  with  respect  to  y-axis. 

Change  in  direction  cosine  of  shock 
normal  with  respect  to  z-axis 

Average  parametric  angle  along  a 
bicharacte  ristic 

Length  of  bicharacteristic  (J) 

Step  size  along  x-axis  between 
successive  data  planes 

Partial  derivatives  of  x,  y  and  z 
coordinates  with  respect  to  the 
N-axis  respectively. 


Coeffi  ients  used  to  compute 

and  iLL 

<9*  dy-  &&  ay-  <9  x 

from  simultaneous  solution  of 
equations  of  the  form 


6~a,  = 


ae 


a* 


L  I  -J 


(see  Subroutine  SETUP) 
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Mnemonic  Variation 


Symbol 


Comme  nts 


ETLAST 

EPLAST 

F(J) 

fnorm 

FB  T( 3,  3) 

FB(3) 

FUNCT(J) 

IADD 

IFIT 

IDP 

NDP 

NFUNCT 

N  OPT 

N  OPTS 
N  O  FN 


3. 


Basis  functions  (J)  used  in  inter¬ 
polating  polynomial 


normalizing  factor  for  unit  normal 
to  body  surface  at  a  point  on  body 

not  used  for  frozen  solution 


not  used  for  frozen  solution 
Solution  for  function  (J)  obtained 
luTacf  mrlaHOn  U8i"6 

it?nDh=  1  in  "fW  ring  of  field  P°ints 
IS  to  be  added  at  shock  IADD  =  0 

otherwise 


Number  of  rings  used  around  center 
ring  in  polynomial  fit.  IFIT  =  1 
always 


Number  of  rings  used  in  surface 
fitting  procedure.  IDP  =  3  always. 

Number  of  rays  used  in  surface 
fitting  procedure.  NDP  =  3  always. 

Number  of  functions  fit.  NFUNCT  -  R 


Number  of  points  included  in  polv 
nomial  surface  fit.  NOPT  =  9 
always 

N  OPT  +1 


,  o - presently 

“  d  ^nalarge  matrix  contains 
everal  different  functions  in 
sequence  (i.  e.  Al) 
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Mnemonic  Variation 


S  ymboi 


Comments 


NFIT 

Number  of  rays  used  around  center 
ray  in  surface  fitting  procedure 

NN 

Iteration  counter 

NOCHAR 

Number  of  bicharacteristics  used 
in  unit  mesh  procedure 

PSIB(J) 

\ 

Flow  angle  at  base  of  bicharacteristic 
(J). 

PSIBT( J,  K) 

Flow  angle  </>  at  base  of  ^character¬ 
istic  (J)  for  set  of  bicharacteristics 
(K). 

PSIAV 

Average  flow  angle  f  along  a 
bicharacteristic 

PPP 

not  used 

PLAST 

Last  value  of  pressure  obtained  in 
iteration  procedure 

PB(J) 

Pressure  at  base  point  of  bicharacter¬ 
istic  (J). 

PI 

not  used 

PPB 

not  used 

OMEGAB 

*>« 

Polar  angle  of  bicharacteristic  base 
point  with  respect  to  y-axis. 

Q1 

Velocity  at  initial  point 

Q232 

Not  used 

QP 

Not  used 

QQB 

Velocity  at  bicharacteristic  base 
point. 

QS 

Not  used 

QB(J) 

Velocity  at  base  point  of  bicharacter- 
istic(J) . 
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Mnemonic  Variation 


Symbol 


Comments 


RB 

R  0 1  ■) 
R02  y 
R03  J 


R1AV  ■) 
R2AV  i 
R3AV  J 

RT1  'J 
RT2  J- 
RT3  J 

RR 

SSB(J) 

SST3 

SSA1 

SS 

SCALE(J) 
THETAB(J) 
THBT(J,  K) 

TA) 

THETAV 


/t0  Radius  frcm  origin  to  bicharacter¬ 

istic  base  point 

Direction  numbers  of  line  of  inter¬ 
section  of  plane  tangent  to  shock 
surface  and  plane  normal  to  shock 
surface  containing  velocity  vector 
behind  the  shock  at  the  initial  shock 
point. 

Average  line  in  shock  surface  for 
locating  new  shock  point  defined 
by  R1AV  =  R01  4  RT1  etc. 

Direction  numbers  of  line  of  inter¬ 
section  of  plane  tangent  to  shock 
surface  and  plane  normal  to  shock 
surface  containing  velocity  vector 
behind  the  shock  at  the  new  shock 
point. 

not  used 

Speed  of  sound  at  base  point  of 
bicharacteristic  (J) 

Speed  of  sound  at  new  data  point 

Speed  of  sound  at  initial  data  point 

Speed  of  sound 

Scaling  factor  for  polynomial  fit 
to  function  (J). 

6g  Flow  angle,  6  ,  at  bicharacteristic 

(J)  base  point 

Flow  angle,  6  ,  at  base  of  bichar¬ 
acteristic  (J)  for  set  of  bicharacter¬ 
istics  (K). 

Temperature  at  initial  data  point 

6  Average  flow  angle,  6  ,  along  a 

bicharacteristic. 
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Mnemonic  Variation 


Symbol 


TAT3 

UAVE 

u 

VAVE 

nr 

WAVE 

Jjr 

US 

LL(r'i 

VS 

WS 

XB(  J) 

XP{4) 

XPT  | 

XNEXT  J 

YB(J) 

** 

YP(4) 

YPT 

ZB(J) 

ZP(4) 

GAMB(J) 

GAMBT(J.K) 

COMMON  LOOP 
"NITER 


Comments 

Temperature  at  new  data  point 

Average  x-component  of  velocity 

Average  y-component.  of  velocity 

Average  z-component  of  velocity 

x-component  of  velocity  behind  shock 

y-component  of  velocity  behind  shock 

z-component  of  velocity  behind  shock 

x-component  of  base  point  for  bi- 
characteristicfJ,) . 

not  used 

x-coordinate  of  new  point  when  plane 
rotation  is  not  employed 

y-coordinate  of  base  point  for 
bicharacteristic  (J). 

not  used 

not  used 

z-coordinate  of  base  point  for 
bicharacteristic  (J) 

not  used 

Specific  heat  ratio  at  base  point 
bichar  icteristic(J). 

Specific  heat  ratio  at  base  point 
of  bicharacteristic  (J)  for  set  of 
bicharacteristics  (K). 

Maximum  number  of  iterations 
allowed 
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Mnemonic  Variation  Symbol 


Comments 


ERROR  =  CHI 


Iteration  limit 


_  5 

1.0  x  10  usually 


COMMON /TTT  APE 
' NREC4 
NREC3 

COMMON /STIPN 
"CSST 

CO  MMO  N  /  B  DYDAT 

RBDY 

XBDY1 

XBDY  2 

XBDY3 

YBDY1 

YBDY2 

YBDY3 

'SLOPE  l  v\ 

‘SLOPE  2  i 
"SLOPE  3  f 

'"bdytyp 

COMMO  N  /O  UT  PNT  / 
IPOINP 


Number  of  records  on  binary  tape  4. 
Number  of  records  on  binary  tape  3. 


Defined  in  input  section 


Spherical  cap  nose  radius  in  cm  or  ft. 

x-coordinate  of  beginning  of  conical 
afterbody  portion  of  body. 

x-coordinate  of  beginning  of  expan¬ 
sion  side  for  elliptical  afterbody 

x-coordinate  of  beginning  of  com¬ 
pression  side  for  elliptical  afterbody 

y-coordinate  of  beginning  of  conical 
afterbody  portion  of  body 

y-coordinate  of  beginning  of  expan¬ 
sion  side  for  elliptical  afterbody 

y-coordinate  of  beginning  of  com¬ 
pression  side  for  elliptical  afterbody 


Defined  in  input  section 


not  used 
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Mnemonic  Variation 


Symbol 


Comment 


SUBROUTINE  INPUT 


"ALPHA  (I) 

BETA(I) 

CH(I) 

CPP(I) 

CM 

0(1) 

GAMEFF 
GAMMA  (I) 

"OMEGA  (N) 

p(D 

Q(l) 

% 

R(D 

"RHO(I) 

SOUND(I) 

'T(I) 


2-D  flow  angle  for  input  point  (I) 
(radians ) 

Mach  angle  for  input  point  (I) 

Total  enthalpy  for  input  point  (I) 
Specific  heat  ratio  for  input  point  (I) 
Mach  number 
Not  used 

Frozen  specific  heat  ratio 

Frozen  specific  heat  ratio  for  input 
point  (I) 

Azimuthal  angle  for  ray  (N) 
Pressure  for  input  point  (I) 

Velocity  for  input  point  (I) 

Radius  for  input  point  (I) 

Density  for  input  point  (I) 

Speed  of  sound  for  input  point  (I) 
Temperature  for  input  point  (I) 


SUBROUTINE  BODY 

D1  "N  Combination  of  coefficients  of 

/  compatibility  equations  used  to 

(  eliminate  pressure  term  and  obtain 

FI  )  one  equation  in  'wo  unknown  (  6  ,  V  ) 

DELT 


Parametric  angle,  6  ,  for  bicharac¬ 

teristic  in  normal  plane  to  body  at 
new  point 


164 


Mnemonic  Variation 


Symbol 


Comment 


DFDPSI 

ar 

9V 

3F 

DFDTH 

99 

DGDPSI 

9  G 
3  ip 

30 

DGDTH 

96 

DENOM 

D THETA 

A6 

DPSI 

a  r 

FN 

G 


SILAST 

THLAST 

X 

SUBROUTINE  BODYFN 
B 

C1  ^ 

Cl 

C3  , 

C1NEW 
C2NEW 
C3NEW 


Partial  derivative  of  function  FN  with 
respect  to  flow  angle  y 

Partial  derivative  of  function  FN  with 
respect  to  flow  angle  9 

Partial  deviative  of  function  G  with 
respect  to  angle  ^ 

Partial  derivative  of  function  G  with 
respect  to  flow  angle  8 

Denominator  for  Newton-Raphson 
method  in  two  variables  (G  and  FN) 
in  terms  of  two  independent  variables 
(  9  ,  V  )■ 

Change  of  angle  9  • 

Change  of  angle  V'  . 

Function  representing  condition  that 
velocity  vector  is  tangent  to  body 

Function  representing  combination  of 
compatibility  equations  eliminating 
pressure,  -fv. 

Previous  value  of  flow  angle  V  in 
iteration 

Previous  value  of  flow  angle  0  in 
ite  ration 

X  =  COS  tf 


Body  function 


Direction  numbers  of  outward  normal 
to  body  in  body  oriented  coordinates 


Direction  numbers  of  outward  normal 
to  body  in  external  coordinate  system 
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Mnemonic  Variation 


Symbol 


Comment 


CALF 

SALF 

FFNORM 


/C&&  (oc)  Cosine  of  angle  between  body  axis 

and  present  external  x-axis 

sAon;  (oc)  Sine  of  angle  between  body  axis  and 

present  external  x-axis 

Norm  of  direction  numbers 


SMAJ 


SMIN 


XNEW 

YNEW 

ZNEW 


Length  of  semi -major  axis  of  body 
ellipse 

Length  of  semi-minor  axis  of  body 
ellipse 


Coordinates  of  body  point  with  respect 
to  body  oriented  Cartesian  system 


Coordinates  of  body  point  with  respect 
to  external  coordinate  system 


SUBROUTINE  CDELTA 

CFSK1 
CFSK2 
CFSK3 

D 1 
D2 
D3 
D4 
D5 
D6 

D7  J 

d  f 

DFNDY  -r1- 

DFNDZ  4^- 


Coefficients  for  quadratic  fit  to  shock 
curve  in  x  =  const,  plane 


Coefficients  of  Mach  forecone  from 
new  shock  point 


Partial  derivative  of  function  FN 
with  respect  to  y-coordinatc 

Partial  derivative  of  function  FN 
with  respect  to  z-coordinate 
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Mnemonic  Variation  Symbol 


Comment 


DGDY 

db 

dt 

Paitial  derivative  of  function  G  with 
respect  to  y-coordinate 

DGDZ 

d& 

di 

Partial  derivative  of  function  G  with 
respect  to  z-coordinate 

DY 

* T 

Change  in  y-coordinate  along 
bicharacteristics 

DZ 

d  z 

Change  in  z-coordinate  along 
bicharacteristics 

DDY 

Change  in  y-coordinate  of  bicharac¬ 
teristics  base  point  as  obtained 
from  Newton- Raphson  procedure  to 
converge  the  base  point  location 

DDZ 

Change  in  z-coordinate  of  bicharac¬ 
teristic  base  point  as  obtained  from 
Newton-Raphson  procedure  to 
converge  the  base  point  location 

FN 

Function  representing  equation  of 

Mach  forecone 

G 

Function  representing  difference  in 
fitted  shock  radius  and  bicharacteristic 
base  point  radius  for  Newton-Raphson 
procedure 

IT  =  IP 

New  shock  point  storage  location 

QAVE 

Average  velocity  along  bicharacteristic 

SSAVE 

Average  sound  speed  along  bicharac¬ 
teristic 

YS(N) 

>  -coordinate  of  shock  point  N. 

ZS(N) 

y  -coordinate  of  shock  point  N 

SUBROUTINE  DIVIDE 

IIPT 

Storage  index  of  every  other  ring  of 

field  points  off  body 
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Mnemonic  Variation 


Symbol 


Comment 


IPTMIN 

SUBROUTINE  FIELD 

C1AV 
C2AV 
C3A  V 

DELXI 

IT 

IP  =  IK 

NP  -  IN 

SUBROUTINE  INTPLT 

RFIT  ) 

OMEG  / 

SUBROUTINE  SECTB 
BT 


CPI  'J 
CP2  > 
CP3  J 

DELY 

DELZ 


Number  of  storage  data  rings 
remaining  after  every  other  ring 
is  dropped 


Average  direction  numbers  for 
streamline  direction 


Length  of  a  birha racte ristic 

Storage  index  for  final  field  point  ring 

Storage  index  for  initial  field  point 
ring 

Storage  index  for  field  point  ray 


Scaled  polar  coordinates  (/u,  ^  )  of 
bicharacteristic  base  point  to  be  used 
in  evaluating  interpolating  polynomial 


Function  representing  difference 
between  actual  body  radius  and  body 
radius  computed  from  present  approx¬ 
imation  to  body  point 

Direction  numbers  of  normal  to 
normal  plane  to  body  containing 
streamline  direction  from  initial 
body  point 

Change  in  y-coordinat.e  of  new  body 
point 

Change  in  ?. -coordinate  of  new  body 
point 


168 


Mnemonic  Variation 


Symbol 


Comment 


FN 

IK 

IP  =  IT 


NP 

SUBROUTINE  SETUP 
AM(J,  1) 

AM(  J,  2) 

AM(J,  3) 

DX' 

DY  - 
DZ 

✓ 

AM(J,  4) 

DELAV 

SDELAV 

C4 

C  5 

C6 

Cl 

C8 

C9 

CIO 

D1  -D9 
SI  BAR ' 

► 

S2BAR 


Function  representing  normal  plane 
to  body  containing  velocity  vector 

Storage  index  for  initial  body  ring 

Storage  index  for  final  body  ring 

(IP=1) 

Storage  index  for  ray  on  which  body 
point  is  located 


Coefficient  of  pressure  in  the 
compatibility  equation  for  bicharac¬ 
teristics  (J) 

Coefficient  of  flow  angle,  6  ,  in  the 
compatibility  equation  for  bicharac¬ 
teristic  (J) 

Coefficient  of  flow  angle,  ^  ,  in  the 
compatibility  equation  for  bicharac¬ 
teristic  (J) 

Distance  matrices  for  coordinate 
directions  along  bicharacteristics 


Coefficient  for  constant  term  in  the 
compatibility  equation  for  bicharac¬ 
teristic  (J) 

Average  parametric  angle  cf 
corresponding  to  bicharacteristic 
curves. 

Sine  of  DELAY 


Average  coefficients  for 
compatibility  equations 


Minors  required  to  solve  for 

de 

-jzg-  ,  etc.  in  EQN  86 
Temporary  storage 
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0 


Mnemonic  Variation 


Symbol 


Comment 


SUBROUTINE  SHOCK 
ALF  i  2 

ALFT3 

ALFT21 

ALFT31 

ALFT22 

ALFT32 

ERRT 

ERROR 


ERROR 1 


IP  -  IK 
IT 

NP 


Present  values  of  direction  cosines 
of  normal  to  shock  with  respect  to 
y  and  z  axis  respectively 

Values  of  ALFT2  and  ALFT3  from 
previous  iteration 


Values  of  ALFT21  and  ALFT31  from 
previous  iteration 

Temporary  storage  for  ERROR  as 
interpreted  in  COMMN/LOOP 

Difference  between  average  pressure 
computed  from  compatibility 
equations  and  the  pressure  computed 
from  shock  relations 

Previous  value  of  ERROR  from  last 
itc  ration 

Storage  index  for  initial  shock  ring 

Storage  index  for  final  shock  point 
ring 

Storage  index  for  ray  on  which  shock 
point  is  located 


SUBROUTINE  STREAM 
DTS 

IK 

IT 

IN 

Q202 

ROINAV 

T  LAST 


Change  in  temperature  along  stream¬ 
line 

Storage  index  for  initial  point  data  ring 

Storage  index  for  final  point  ring 

Sic  rage  index  for  ray  on  which  field 
point  or  body  point  is  located 

Velocity  squared  divided  by  2 

Inverse  of  average  density  along 
streamline 

Initial  value  of  temperature  along 
streamline 
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Mnemonic  Variation 


Symbol 


Comment 


SUBROUTINE  STOP 
NSTOP 


NFILE 


Integer  test  constant  for  stopping 
integration  if  X  >XMAX  or 
NPLANE  >■  NPLMAX 

Number  of  data  planes  stored  on 
binary  tape  2 


SUBROUTINE  FLOW 

CPAV 

DHDP 

DUN 

ERR 

HT 

HPREV 

PREVP 

QS  COMMON  AFT3D 

TIP 

UNAM 

UNIP 

VTIP 


Frozen  specific  heat  behind  shock 

Change  in  total  enthalpy  divided  by 
change  in  pressure 

Change  in  velocity  component  normal 
to  shock,  across  shock 

Difference  between  free  stream, 
total  enthalpy  and  computed  total 
enthalpy  across  shock 

Total  enthalpy  behind  shock 

Previous  value  of  HT  for  previous 
ite  ration 

Previous  value  of  pressure  behind 
shock  from  previous  iteration 

Velocity  behind  shock 

Temperature  behind  shock 

Free  stream  velocity  component 
normal  to  shock 


Component  of  velocity  behind  shock 
in  direction  normal  to  shock 

Component  of  velocity  behind  shock 
in  direction  tangent  to  shock 


x,  y  and  z  velocity  components 
upstream  of  shock 


USP 

VSP 

WSP 
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Mnemonic  Variation  Symbol 


Comment 


US  , 

I  y  and  z  velocity  components 

VS  COMMON  AFT3D  across  shock 

WS) 


1  72 


A.  7  LOGICAL  FLOW  DIAGRAMS 

Detailed  flow  diagrams  including  all  the  equations  and  symbols 
are  probably  not  so  valuable  as  logical  charts  when  the  program  is  written 
in  FORTRAN  IV  as  is  the  case  here.  For  this  reason  only  the  general 
logical  flow  diagrams  are  presented  with  no  attempt  being  made  to  be 
complete  in  every  detail. 
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Appendix  B 

CALCULATIONS  OF  FIELD,  BODY,  AND  SHOCK  POINTS 
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TABLE  B-I 


CALCULATIONS  AT  A  FIELD  POINT 
Step  Size  Ax  =  0.  01460736 


Initial  Point  Data:  (All  Quantities  Nondimensionalized) 


JL  JL 

-0.45340951  0.75886150 


JL 

0. 21346586 


0.  62946524 


9  (Radians) 
0.  3839265 


Final  Point 


First 

Iteration 

Second 

Iteration 

Third 

Iteration 

x: 

-0.  43880215 

-0.  43880215 

-0. 43880215 

0.  76503433 

0.  7649  3999 

0.  76493883 

r- 

0. 55582996 

0. 55576114 

0. 55576062 

0.  20518689 

0. 20505461 

0. 20505998 

r- 

0. 63728678 

0.  63545944 

0. 63556746 

9: 

0.  37413336 

0. 37405871 

0. 37406140 

V- 

0. 28928342 

0. 28921835 

0. 2892211 3 

Convergence  Test  Criterion 


(r)  (tt-0 

f  -_P 


0. 00001 


0.  55134514 

JL  (Radians) 
0.  2978896 

Fourth 

Iteration 

-0.43880215 
0.  76493885 
0.  55576064 
0.  205059835 
0. 635562873 
0. 37406141 
0.  28922114 
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I 


TABLE  B-II 


CALCULATIONS  AT  A  BODY  POINT 
Step  Size  Ax  =  0.  01460736 


Initial 

Point  Data: 

X 

-0. 45340951 

y 

-0.  44568498 

y 

0.  77194911 

jP_ 

0.  19598220 

±_  e_ 

0.61962985  -0.  2285j  .11 

V 

0.  41459165 

Final  Body  Point: 

First 

Iteration 

Second 

Iteration 

Third 

Iteration 

X 

-0. 043880215 

-0.  043880215 

-0.  043880215 

y 

-0.  44928833 

-0.  44929017 

-0.  4492910 

9 

0. 77819028 

0.  77819325 

0. 77819480 

f> 

0.  1868241 1 

0.  18683943 

0. 18683876 

j 

9 

0. 62771747 

0.  62712205 

0. 627113444 

1 

9 

-0. 2212261  3 

-0. 22119960 

-0.  221  199  38 

f 

0. 40008304 

0.  40009536 

0. 40009460 
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TABLE  B-III 


CALCULATIONS  AT  A  SHOCK  POINT 
Step  Size  Ax  =  0.  01460736 


Initial  Point  Data: 

& 

-0.  45340951 


_y  ± 

1.2645526  0.0 


_P 

0.  34633269 


± 

0.  77348271 


e_ 

0.53874810  0.0 


-0.  64625884 


76311828 


0.  0 


where  N  -  normal  vector 
to  shock 


Final  Point: 


First 

Second 

Third 

Iteration 

Iteration 

Iteration 

X 

-0. 43880215 

-0. 43880215 

-0.  43880215 

V 

1.  2769231 

1. 2764755 

1. 2768739 

0.  0 

0.  0 

0.  0 

0.  1  6793336 

0.  1721  8783 

0.  1  5803909 

fs 

0.  14470096 

0. 1 3455693 

0.  16085586 

9 

0.  21948106 

0.  22696227 

0.  20741885 

9 

0.  52119923 

0.  511  89  300 

0.  53575372 

0.  0 

0.  0 

0.  0 

-0.  61823331 

-0. 64325045 

-0. 64156076 

N, 

0.  78599464 

0. 76565583 

0. 76707221 

N> 

0.  0 

0.  0 

0.  0 
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TABLE  B-III  (Cont. ) 


Forth  Iteration 

X 

-0.43880215 

y 

1. 2768463 

} 

0.  0 

Pc 

0.  15898510 

Ps 

0.  15904675 

? 

0. 20877892 

e 

0. 53413956 

V 

0.  0 

A '* 

-0.  6-H  52290 

Ny 

0.  76710389 

N> 

Convergence  Test  Criterion 


(ft)  (ft) 


Pc 


Ps 


Ps 

(n) 


< 


0. 00001 


Fifth  Iteration 
-0. 43880215 
1. 2768457 
0.  0 

0.  15900633 
0. 15900626 
0. 20880932 
0. 53410338 
0.  0 

-0. 64152295 
0. 76710384 
0.  0 


where 

Pc  =  average  pressure  from  compatibility  equation 
Ps  -  pressure  from  shock  relations 
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Appendix  C 

AUTOMATED  INPUT  PREPARATION 


For  certain  configurations,  the  flow  over  the  nose  region  may  be 
described  by  two  space  variables.  In  these  cases,  the  preparation  of  input 
data  for  the  present  three-dimensional  program  will  be  facilitated  by 
certain  modifications  of  the  axisymmetric  program  described  in  Reference 
60.  For  convenience,  FORTRAN  listings  of  the  relevant  portions  of  the 
referenced  program  are  reproduced.  The  shaded  statements  indicate  the 
changes. 
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BEST  AVAILABLE  COPY 
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The  axially  symmetric  program  of  Reference  60  ordinarily  uses 
a  mesh  composed  of  streamlines  and  rays  normal  to  the  body  surface.  In 
order  to  prepare  input  data  for  the  present  program  it  is  desirable  to  rotate 
these  rays  about  a  fixed  point  on  the  body  surface  until  they  are  normal  to 
the  wind  axis.  In  doing  so,  however,  it  is  essential  that  the  resulting  ray 
be  located  so  that  the  backward  characteristics  from  the  first  plane 
computed  by  the  3-D  program  intersect  the  initial  data  plane.  This 
condition  is  equivalent  to  requiring  that  the  angle 

SI  +  0  <T  90° 

The  indicated  program  changes  will  rotate  the  rays  through 
successive  1 -degree  increments  as  soon  as  it  is  safe  to  do  so  (  SI  <  88.  5°) 
and  will  then  punch  the  necessary  data  in  the  proper  format.  FORTRAN 
unit  9  is  used  in  the  program  as  written  above,  but  this  may,  of  course, 
be  altered  to  suit  the  requirements  of  the  local  computing  system. 

C.  2  INPUT  PROCEED! RE 

The  input  to  the  axially  symmetric  program  is  set  up  as  described 
in  Reference  60,  VOL.  II.  If  the  ray  rotation  and  punching  are  desired,  the 
user  must  punch  a  1  in  column  20  of  "Card  Y"  (Ref  60,  p.  14  Vol.  II  p.  14). 
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